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Abstract 


Software was developed to construct NURBS curves that approximate the geometries of iced airfoils. 
Users specify a tolerance that determines the extent to which the approximating curve follows the rough 
ice. This ability to smooth the ice geometry in a controlled manner will assist the generation of grids 
suitable for numerical aerodynamic simulations, and ultimately aid studies of the effects of smoothing 
upon the aerodynamics of iced airfoils. The software was applied to several different types of iced airfoil 
data collected in the Icing Research Tunnel at NASA Glenn Research Center, and was found to efficiently 
generate suitable approximating NURBS curves for all geometries. This method is an improvement over 
the current “control point formulation” of Smagglce (v. 1.2). In this report, we present the relevant theory 
of approximating NURBS curves and discuss typical results of the software. 


Introduction 

The detrimental effect of ice upon airfoil performance has long been a safety concern, and is largely 
studied using experimental models in special wind tunnels such as the National Aeronautics and Space 
Administration’s Icing Research Tunnel at the Glenn Research Center. In recent years, the reduced 
aerodynamic performance of airfoils with moderate ice has been successfully simulated via computational 
fluid dynamics (CFD) [1-4]. There is hope that CFD can be used to simulate aerodynamic performance 
under more general icing conditions. Experimental and computational efforts would then serve more 
complementary roles in the development of safe aircraft, and safety margins could be enhanced at reduced 
costs. 

In general, the presence of ice on an airfoil greatly increases the difficulty of a CFD analysis over that of a 
clean airfoil. The ice often has deep narrow crevices and exhibits varying degrees of roughness. The 
associated flow over iced airfoils is very complex. In particular, existing grid generation codes are 
generally ill equipped to accommodate the wide variety of shapes and sizes of ice that is typically found 
on an airfoil. So instead of being automated as it is for a clean airfoil, the generation of quality grids for 
iced airfoils is frequently a labor intensive, interactive process. To address this problem, NASA has been 
developing an interactive software toolkit [5] called Smagglce 2D. This software enables the user in the 
tasks of geometry preparation, domain decomposition, block boundary discretization, gridding, and 
linking to the flow solver for a two-dimensional airfoil. 

The current version (vl.2) of Smagglce 2D provides an option to represent the ice geometry via a Non- 
Uniform Rational B-Spline (NURBS) curve [6-8]. NURBS have a number of attractive features, 
including fast and numerically stable algorithms, and easy-to-understand geometric interpretations. 
Designers of the current software found it convenient to place the control points of the NURBS curve on 
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the x-y data that specifies the ice geometry. Ordinarily, control points are not placed on the geometry to 
be represented, but in this case the error is small due to the large number and high density of the given 
data. Nodes are later placed along the representation for numerical simulation. 

Smoothing of the ice geometry representation is often required to generate a quality grid by state-of-the- 
art grid generators. Smagglce 2D permits the user to smooth the NURBS representation by reducing the 
number of control points [5]. In this method, the user can only coarsely control the smoothed NURBS 
representation of the ice geometry. This smoothing facilitates the subsequent gridding and numerical 
simulation, but the lack of control makes it difficult to access the effects of smoothing on the CFD results. 
In our approach, the user requires the NURBS curve to satisfy a tolerance, i.e., a characteristic distance 
between the ice geometry and its representation. In this way, a user can more precisely control the extent 
to which the representation is smoothed. Not only will this ease the difficulty of generating quality grids 
over ice, but it will also aid investigation of the relation between ice roughness and aerodynamic 
characteristics of the iced airfoil. 

We developed a software package that uses a NURBS curve to represent a two-dimensional cross-section 
of the ice surface as closely as desired by specification of a tolerance. Given a set of point data that 
defines the two-dimensional ice geometry, the software finds a NURBS representation to within the 
specified tolerance, which may be expressed either as a maximum distance or as a maximum root-mean- 
square (rms) distance between equivalent points of the curve and the prescribed data. The user is given 
various suitable options that further modify the NURBS representation, including the degree of the 
NURBS basis functions and whether or not to fix end-point derivatives. The latter option was found to 
permit the achievement of tighter tolerances. The code was applied to a wide range of experimental data 
from the Icing Research Tunnel, and was found both robust and efficient in all cases. In this report, we 
further document this software by providing requisite theory and typical results. The FORTRAN 77 code, 
which is included as an appendix, is available for inclusion in the next version of Smagglce 2D (vl.8). 

Theory 

We here review the theory of Non-Uniform-Rational-B-Splines. Our emphasis is upon generating a curve 
C (u) = (x(u),y(u)) that approximates prescribed data points Q A (k = 0,1,..., mdata) in two dimensions. 
Both the theory presented below and the software are limited to a large subclass of NURBS curves, the 
non-rational B-spline curves, to avoid nonlinearities associated with determining all parameters needed in 
the more general class. B-spline curves (we henceforth drop the word non-rational) cannot exactly 
represent certain curves (e.g., perfect circles) that a more general NURBS can. Nevertheless, a B-spline 
can represent any smooth curve to within a tolerance, which is all that is needed in the present application. 
A list of symbols is provided in Appendix A. 

Preliminaries 

We represent a B-spline curve as the finite sum, 

C(u) = f j N. p (u)P i , (1) 

i=0 

with curve parameter u in the interval [0,1] , P the control points in a two-dimensional space, and 
A,, (n) the /;th degree B-spline basis functions. These basis functions are defined on a knot vector 

U — {u 0 ,..., u mknot } , 

where u j (i = ()...., mknot) are the knots. For present puiposes, the knot vector has the form 
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( 2 ) 


U = j Q.o, ....o , U p+l ,u p+2 ,..., u mknol _ p _t , IUJ ^ , 

[ p+i p+ 1 J 

in which the unknown knots increase in numerical order in the open interval (0,1). Knowledge of the 
degree p (which the user provides), the knot vector U, and the control points P is required to fully 
specify a B-spline. In the procedure to be described below, a sequence of approximating B-splines is 
generated. When the procedure is successful, each spline in the sequence more closely represents the 
given data. 

Piegl and Tiller [6] furnish several properties of the basis functions and B-spline curves. We here list 
those that are most relevant to the present application, and refer the interested reader to Piegl and Tiller 
for additional details and references. 

Non-zero basis functions: According to the local support property, Ni, p (u ) = 0 for u outside the interval 
[ u, , iii+p+i ). It follows that for any given value of the curve parameter u, at most only p + 1 of the n + I 
basis functions that appear in Equation (1) are non-zero. Therefore, for any value of u, only p + I 
consecutive terms in the equation actually need to be determined and summed. 

Efficient computation of basis functions: Given a value u of the curve parameter and a knot vector U, 
the non-zero basis functions may be computed simultaneously and efficiently. 

Relation between the number of knots, the number of terms in the sum, and the degree: These three 
quantities are related by the simple equation 

mknot = n + p + 1, (3) 


after subtraction of one from both sides. In the iterative procedure discussed below, the degree p is fixed, 
and mknot is increased so that the B-spline may more closely follow the ice geometry. This relation 
requires an identical increase in n. 

Continuity and differentiability of C(u): Given a knot vector of the assumed form, the associated B- 
spline curve is continuous, infinitely differentiable in the interior of knot intervals, and at least p times 
differentiable at each knot. The curvature of the B-spline is often required when determining the 
placement of nodes along the curve for grid generation [9]. Calculation of curvature involves the second 
derivative C"(u). Continuity of curvature is assured if p> 3. 

Endpoint interpolation of C (w): At a = 0 all basis functions are identically zero except N (>.,,(()) = 1; 
similarly, at u = 1 all basis functions are identically zero except /V„ ; , (0) = 1. It follows from Equation (1) 
that a B-spline must interpolate the first and last control points: C(0) = P 0 and C(l) = P n . All B-splines of 
interest to us also interpolate the first and last points of the prescribed ice geometry. This immediately 
leads to the conclusion that for our purposes first and last control points respectfully coincide with the 
first and last data points: 


p 0 =Q 
p„=Q 


0 

mdata * 


(4) 


Endpoint derivatives of C (u): First-order derivatives C'(0) and C'( 1 ) of a B-spline are related to the 
first and last pairs of control points via the respective expressions [6] 
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( 5 ) 


C'(0) = — ^-(P,-P 0 ) 

U P+ 1 

C’(D = — — 

M'mknot- p - 1 

These equations may be solved for control points P l and P ;1 : 

Pi =Q 0 +— C'(0) 
p 

P »-,=Q^-— C'(D 

p 

The above relations will be useful in the development of B -spline representations that satisfy prescribed 
endpoint derivative conditions. 


Global Approximation by a B-Spline 

Piegl and Tiller [6] offer a procedure for the global approximation of prescribed data by a B-spline curve 
provided the endpoint derivatives are free. We summarize their method below, and later extend it to the 
case for which endpoint derivatives are specified. In both cases, a suitable global approximation is 
obtained by iteration. The user provides certain information such as the point data to be approximated, a 
tolerance to be satisfied, and an initial value of n. Based upon this information, an initial approximating 
B-spline is developed using the method of least-squares. This approximation is tested against the 
tolerance requirement. If the requirement is met, the approximation is accepted. If not, the number of 
knots and terms in Equation (1) is increased and a new B-spline is generated. This process continues until 
the specified tolerance is achieved or the process fails. Failure could occur, e.g., if the specified tolerance 
is too small. 

Free Endpoint Derivatives 

The data Q k , degree p of the approximating B-spline, an initial value of n, and a tolerance requirement 
are presumed given. Because we seek an approximating B-spline that interpolates first and last data 
points, first and last control points are given by Equation (4). The remaining control points and the knot 
vector are obtained as follows. We first parameterize the given data; i.e., a parameter value is assigned to 
each data point. While there are several ways in which this can be done, the chord length method is 
generally satisfactory. Let d be the total chord length, 

m data — 1 

d= £ |Q t+1 -Q*|, (7) 

k = o 


and require first and last points to correspond to the respective parameter values u 0 = 0 and u mdata = 1 . 
The remaining data parameters are then defined by the recursive equation 



for k = 1, . . . , mdata - 1 . These data parameters remain unchanged throughout the rest of the analysis. Data 
points can be numbered from either end of the ice geometry. To be definite, we assume that points are 
numbered consecutively from the upper to lower airfoil surface. 
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An initial knot vector must be generated that in some sense corresponds to the distribution of data 
parameters. Recall from Equation (2) that the first and last p + 1 knots are identically zero and unity, 
respectively. To determine the remaining knots, we define a new temporary variable d via the expression 


mdata + 1 
n - p + 1 


(9) 


which represents the average number of data points per knot span of nonzero length. The remaining knots 
are found from the following equations [6] : 

i = int (j * d) 

a = j*d-i (10) 

u p+j = (l-a)*u M +a*u i for j = 1 p 

Here, the staircase function y = int(v) is the largest integer such that y < x. The above definition is 
successful in placing at least one data parameter u k in every knot span of nonzero length if d > 1. 
Provided this is true, a key symmetric matrix ( N 1 N. which is encountered below) in the approximation 
procedure is positive definite and well-conditioned [6-8], at least in theory. 

Determination of the internal control points P ( (i = I, - I ) is the final major step in defining the initial 
B-spline. These control points are selected so as to minimize the sum of the squared distances of the curve 
from data points Q,,..., Q : 

m data - 1 

I |Qr-C07,)f 

k = 1 


Take the derivative of this expression with respect to control point P ( and set the result to zero. This 
operation yields n — 1 equations that can be collectively represented by the matrix equation 

(N t N)P = R. (11) 


Here, N is the ( mdata - 1) x (n - 1) matrix 


N = 


v ^ ] , p mdata— 1 ) ^ n—\,p mdata— 1 i 


( 12 ) 


J 


whose elements consist of basis functions evaluated at certain values of the data parameters. Recall that 
no more than p + 1 of the basis functions are non-zero for any given value of the curve parameter u. 
Thus, each row of N contains at most p + 1 non-zero entries. Also appealing in Equation (1 1) is a 
(n - 1) x 2 matrix P of unknown control points, 


P = 


V P «-1 J 


(13) 


and a (n - 1) x 2 matrix R: 
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f mdata—l 


R = 


Z N hp( u k) R k 


Z N n - !.,(«* 

V *=i 


(14) 


Each (1x2) row vector R A . that appears above is defined by the expression 

R* =Q, -N 0 , p (u k )Q 0 -N n Ju k )Q mdata (15) 

with k = 1 , . . . , mdata - 1 . (Recall that Q k = ( x k ,y k ) is defined within a two-dimensional space, which 
explains why P and R each have two columns and why R A is a two-element vector.) Equation (11) thus 
represents two 1 i near systems of equations having the same matrix coefficient but with different right- 
hand sides. 


Given that each knot span of non-zero length contains at least one data parameter, the (n - 1) x (n - 1) 
matrix N r N (with superscript T denoting the transpose operation) is symmetric, positive definite, and in 
principle well-conditioned. (However, in practice the matrix can become singular as n becomes large and 
d in Equation (9) approaches unity. This is presumably due to limitations of finite arithmetic on a digital 
computer.) Moreover, it has a semi-bandwidth of p + 1 , and can be stored in a compact form. The 
Cholesky method efficiently factorizes the matrix. Control points P, , . . . , P n l are then easily found via 
back substitution. 


Solution of Equation (11) leads to a complete B-spline representation of the ice data. The next major step 
is to determine whether or not this representation satisfies the specified tolerance requirement. Two such 
tolerance criteria are implemented. The user may either specify a maximum tolerance requirement, 

<ax = max |C (u k ) - Q ( . | < £ max , (16) 

\<k<mdata -\ ' 


or a rms tolerance requirement: 


d 


rms 


m data — 1 

Z IQ, -CM 2 

— kd < 

mdata—l 


(17) 


In the above, c/ m . ix is the maximum separation distance of all the distances |C(tq) — Q, | , £ max the 
maximum tolerance, r/ rms the rms distance, and £ ms the rms tolerance. If the specified tolerance 
requirement is satisfied, the B-spline is accepted as a suitable approximation to the ice geometry. If not, a 
new knot vector having additional knots is created, and a new B-spline approximation is generated using 
the above procedure. This cycle is repeated until the desired tolerance requirement is satisfied. Details 
concerning how the new knot vector is generated and how the tolerance criteria are implemented are 
reserved for a later section. 


Before leaving this section, we note that C(u k ) is generally not the closest point on the curve to data 
point Q k . Nevertheless, the distance COT, ) - Q, , which appears in both tolerance criteria, is useful and 
convenient for the intended use. An alternative, more exact maximum tolerance criterion, 


max 

\<k<mdata—\ 


C (K ),n,n)-Q, ^ C, 


could be used. Here u = (u k ) mm is the parameter value that minimizes the distance |G(w) - Q k | for the kth 
point. Determination of each u = (u k ) mm is a nonlinear problem that can be solved, e.g., using a Newton 
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iteration procedure [6]. The advantage in using the latter tolerance requirement in the present application 
does not appeal - to offset the extra required computational effort: up to mdata - 1 values of u - (u k ) mm 
are required per loop of the iteration that ultimately results in an acceptable B-spline curve. In any case, 
because |C((iq ) mm ) -Q A | < C(/7, ) - Q t |, we are assured that any curve that satisfies Equation (16) also 
satisfies the above exact requirement. Our decision to use u = u k rather than u = (u k ) mm may result in a 
larger final value of n than that needed to meet the above exact tolerance requirement. 

Fixed Endpoint Derivatives 


Fixing of the endpoint derivatives requires a straightforward modification of the above least-squares 
procedure. Recall from Equation (5) that the endpoint derivatives of a B-spline curve are related to the 
first and last pairs of control points. If we specify these derivatives, Equation (6) gives control points 
Pj and P n _j so that P () , P, , P^, , and P n are all known prior minimizing the sum of squared distances. In this 
case, we minimize the sum with respect to the remaining control points: P,,...,P n _ 2 . The least-squares 
procedure is altered only in that matrices N, P, R, and vector R A must be redefined. If endpoint 
derivatives are fixed, N is the ( mdata — 1) X (n — 3) matrix 


^2,p( U d ••• N n -2' P ( U l) 


N = 


X mdata— 1 ) ^ n-2, p (F mdata- 1 ) y 


(18) 


and {n - 3) x 2 matrices P and R have the respective new definitions 


P = 


f p A 

r 2 


\ P n-2y 


(19) 


and 


f 


R = 


m data — 1 ^ 

£ N 2, p (M k ) R k 


k = 1 


m data - 1 

£ N n-2,p(K) P k 


V A=1 


J 


( 20 ) 


Finally, each vector R k is defined by 

R*=Q*-( N o,p « ) p 0 + N 1 , P (u k )P| + N n _ lp (u k )P_, + N np (u k )P„ ) , 


( 21 ) 


with P 0 ,P p P n j, and P n given by Equations (4) and (6). The unknown control points are found by 
solving matrix Equation (11) with these modified definitions. 


Software Implementation 

In this section, we discuss how the two tolerance requirements are implemented in the software, and 
provide formulas for calculation of default values of endpoint derivatives. This background will prove 
important in understanding sample results produced by the code as discussed in the next major section. 
Brief descriptions of major elements of the code may be found in Appendix B, and the code itself is given 
in Appendix C. 
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Maximum Tolerance 

The maximum tolerance requirement, Equation (16), is implemented as follows: Beginning with k = 1, 
the distance |C {u k ) - Q t | is compared with the maximum tolerance. If the distance is less than or equal 
to £ max , k is incremented by one and the test is repeated. If the test fails for any k, the knot span associated 
with u k (i.e., the span such that u i < u k < u i+l ) is labeled as non-converging. To avoid unnecessary 
evaluations, the remaining distances within that knot span are skipped. The next distance compared with 
^max i s that corresponding to the first u, located within the next knot span. The test continues until all 
knot spans are labeled as converging or non-converging. 

If the specified tolerance is not achieved by all knot spans, a new knot vector is generated as follows: A 
new knot is inserted at the midpoint of all non-converging knot spans, and both n and mknot are increased 
by the number of added knots. Next, the knot vector is inspected to insure that each knot span of nonzero 
length contains at least one data parameter u k . If it does, the above least-squares procedure is repeated to 
generate a new set of internal control points and a new B-spline curve. The maximum tolerance criterion 
above is checked again. If the maximum distance <7 max from the prescribed data is less than or equal to the 
specified tolerance £ max , the iterative procedure ends. If the tolerance is not achieved, a new knot vector 
is generated by inserting a new knot at the midpoint of each new non-converging knot span, and the cycle 
repeats. 

A knot vector formed by repeated insertion of knots may have one or more knot spans of nonzero length 
that are empty, i.e., ones that do not contain at least one data point parameter u k . In this case, the matrix 
N 1 N is no longer guaranteed to be positive definite and well-conditioned. Provided the new value of n 
is sufficiently small ( n < mdcitci + p -1 ), a completely new knot vector is generated using Equation (10) 
that has at least one data point parameter located within each knot span of nonzero length. A new B-spline 
is generated based upon this new knot vector, and the iteration procedure continues as usual. If the new 
value of n is too large ( n > mdata + p - 1 ), the procedure will abort. 

Root-Mean-Square Tolerance 

The maximum tolerance requirement of Equation (16) simply requires that the maximum distance that 
separates the ice geometry from equivalent points on the B-spline curve is less than a specified distance 
called the maximum tolerance. This measure of closeness of the spline to the ice geometry does not 
distinguish between a curve that is close almost everywhere with the exception of one or a few values of 
u k and one that is far (but less than or equal to the maximum tolerance) almost everywhere. The ability to 
specify a rms tolerance, £' rms , permits the user to require the curve to be close to the data in a rms sense. 
We also point out that the rms distanced^, which is defined in Equation (17), is closely related to the 
actual quantity that is minimized in the least-squares procedure. 

The iterative procedure described above to achieve the maximum tolerance criteria is modified to satisfy 
the rms requirement. In this case, the user must supply an initial maximum tolerance £ m . a , a rms 
tolerance £ nTis , and a scalar a. The rms tolerance is restricted to the interval (0, £ max ) while cnnust lie in 
the interval (0,1); the default value given a is 0.9. During iterations, the initial maximum tolerance 
requirement is satisfied as described above. Prior to exiting the loop, the rms tolerance test is performed. 

If <7 rms < £' rms , the loop exits normally. Otherwise, the maximum tolerance is reduced according to the 
expression £ max =a* d max and each knot span is re-examined to determine the non-converging knot spans 
relative to the new maximum tolerance. Knots are inserted at the midpoint of each such span, and 
iterations continue until the new maximum tolerance requirement is satisfied. The rms criterion is then 
tested again. The loop exits normally when both the modified maximum tolerance and the rms tolerance 
requirements are met. 
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Endpoint Derivatives 

First-order endpoint derivatives are specified in polar coordinates. The angle must be measured from the 
positive v axis and specified in units of degrees. Because these derivatives are vectors, both direction and 
magnitude must be supplied for each endpoint. For convenience, default values for these derivatives are 
calculated based upon finite differences of the first and last pairs of data points: 


C'(0) = 


C’(l) = 


Qi -Q 0 _ Qi -Q„ 


Vt i U Q Mj 

Q mdata ^^indata - 1 . 


1 — u„ 


( 22 ) 


Before leaving this section, we note that these derivatives are taken with respect to the curve parameter u, 
and whether this parameter is increasing or decreasing when moving along the B -spline ultimately 
depends upon the sense of direction of the original data. The user must take this into account when 
specifying non-default values for C'(0) and C'(l). 

Application: Sample Results for an Ice Geometry 

The FORTRAN 77 code was compiled and run on a 1.1 GFlz personal computer (Intel Pentium III 
processor under Linux). Suitable B-spline approximations that corresponded to a wide range of tolerances 
were obtained for three ice geometries, which were provided by the Icing Research Tunnel. To illustrate 
typical behavior, we investigate several representations for one ice geometry. 

Most of our runs - and all discussed here - used cubic B-splines. The few runs conducted with degree 
p = 4or 5 yielded B-splines that, on a graph at least, appeared equivalent to a cubic B-spline for the same 
geometry and parametric values. No novel features were discovered in these runs. 

Figure 1 shows the front portion of a clean airfoil, the attached ice, and an approximating B-spline. In this 
case, we required the B-spline to satisfy a maximum tolerance of £ max = 1.0xl0“ 2 . Endpoint derivatives 
were free, and n was set initially equal to the default value of three. The original data was normalized 
with respect to the chord length; all lengths, including tolerances, are therefore similarly normalized in 
this and all subsequent figures. A total of 525 points defines the full ice geometry. A significant fraction 
of these points are denoted in the figure by small plus signs. Though in this figure one cannot determine 
the particular point C(u k ) on the B-spline curve that corresponds to the Arth data point Q t , it appears that 
every point on the B-spline curve is located within the maximum distance l.Ox 1(T 2 of some point on the 
ice geometry. 
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Figure 1: Global view of front portion of clean airfoil, ice data, and approximating 
B-spline {p = 3; £ mm =1.0xlCT 2 ; free endpoint derivatives; initial n = 3) 



Figure 2: Close-up view of three B-spline representations of ice data at the indicated moderate 
values of the maximum tolerance ( p = 3; free endpoint derivatives; initial n = 3) 
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Figure 2 shows a magnified section of the clean airfoil, the ice data, and three approximating splines near 
the upper icing limit. Maximum tolerances of the approximating curves are, in decreasing order, 
f max = I .Ox 10 ’,5. Ox 10 ', and 1.0x10 '. Visual examination shows that the average distance between a 
B-spline and the ice data correlates with the specified maximum tolerance. At this scale one can see that 
all three B-spline curves interpolate the first point Q 0 of the ice data as required. The consequence of not 
specifying the endpoint derivatives is also seen: the slopes of the three splines differ enormously at the 
upper endpoint. 

The beginning of the ice data and a B-spline (the one with the smallest tolerance in Figure 2) is shown at 

greater magnification in Figure 3. Circles on the B-spline curve represent C(u k ) for k = 0 1 0. Observe 

that the curve passes somewhat closer to each point Q A (indicated by plus signs) than the nominal 
separation distance \C(u k ) - Q, | . 



Figure 3: Magnified view near upper icing limit showing relation between 
C(u k ) and Q k ( p = 3; free endpoint derivatives; £ mnx = 0.001; initial n = 3 ) 


The two B-spline curves that appear in Figure 4 were generated with the maximum tolerance set at 
£ mwi = 1 .0x 1 0 The curve with large wiggles near the endpoint was constructed with free endpoint 
derivatives; the other curve had its endpoint derivatives specified as the default values. Both curves were 
generated using the appropriate default initial values for n (n = 3 for the free endpoint derivative curve; 
n = 4 for the fixed endpoint derivative curve). The success of the second curve as a suitable representation 
of the ice geometry is due to specification of the endpoint derivatives. (One might argue that the different 
initial values used for n might also contribute to the observed difference in these B-splines because the 
initial knot vectors are different. Flowever, for this geometry we found that the same B-spline is obtained 
in the free endpoint derivative case whether an initial value of n = 3 or n = 4 is used. This is because at 


NAS A/CR— 2004-2 13071 


11 




0.014 


0.012 

0.01 


0.008 


0.002 0.004 0.006 0.008 0.01 
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Figure 4: Effect of fixing endpoint derivatives in a B-spline 
( p = 3; £ max = 0.0001; initial n = 4 ) at small tolerances 




Figure 5: Close-up view of three B-spline representations of ice data at the indicated moderate 
values of the maximum tolerance ( p = 3; fixed endpoint derivatives; initial n = 3) 
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some stage of the iterative process, both procedures use Equation (10) to form a new knot vector for the 
same value of n. All subsequent steps in the two procedures to generate the final B-splines are then 
identical.) This same behavior was found in corresponding representations for both of the other ice 
geometries investigated. It is widely known that B-splines have a tendency to exhibit wiggles when 
approximating noisy data at small tolerances [6]. We note that, despite specification of the endpoint 
derivatives, unacceptably large wiggles in the B -spline representation may arise at smaller tolerance 
levels. These were observed for the present geometry with a maximum tolerance of £ mm = 2.5 x I (T : . 

It may be advantageous to fix the endpoint derivatives at more moderate levels of tolerance. For example, 
the slopes of the iced airfoil near the icing limits are sometimes relevant in the merging of the ice 
representation with the clean airfoil to avoid the introduction of significant discontinuities [9] . In any 
case, specification of endpoint derivatives can have a dramatic effect on the entire B-spline, not just near 
the endpoints. Figure 5 shows three B-splines for which the endpoint derivatives were set to their default 
values; all curves in the figure appear to have the appropriate endpoint derivative. The maximum 
tolerances of the three splines in Figure 5 are the same as those in Figure 2; corresponding curves in the 
two figures should be compared. These two figures, along with Figure 4, suggest the following two 
general rules: a) the effect of specifying an endpoint derivative decreases with distance (along the curve) 
from the endpoint, and b) the distance over which the endpoint derivative has a significant effect is 
inversely related to the tolerance. If the tolerance is sufficiently large, as it is for the two £ max = I .Ox 1 0 2 
curves, the curves can appear significantly different from each other for the whole domain. Furthermore, 
by controlling the values of the endpoint derivatives, one can obtain intermediate representations of the 
ice geometry. These observations appeal - consistent with the local support property of the basis functions 
and the (frequently) inverse relationship between tolerance and the number of knot spans in the resulting 
B-spline. 

We note that in those runs reported above with acceptable representations of the ice geometry (i.e., no 
significant wiggles), the ratio d nvM / £ max varied over the range (0.77, 0.99). Moreover, the variation of this 
ratio was not monotonic with respect to the tolerance. That the ratio is not constant should not be 
unexpected because the algorithm does not control this quantity beyond the requirement that it be located 
within the interval (0,1]. Similarly, the algorithm does not control the ratio d rms / £ rms when the rms 
tolerance is specified; it too is limited to the interval (0,1], The variation of these ratios has an interesting 
potential consequence when two B-spline approximations of the same ice geometry is generated with 
either the maximum or rms tolerances slightly different but all other parameters unchanged. The usual 
expectation is that the distance d mia or r/ nns , depending upon which tolerance is specified, will be smaller 
for the spline with smaller specified tolerance. While this pattern frequently will be true, there may be 
exceptions (e.g., r/ niax larger for the spline generated with smaller £ mdX ) because the corresponding ratio is 
not constant. This behavior in no way detracts from the anticipated use of the software, which is to easily 
generate in a controlled manner B-spline representations of ice geometries that have a wide variety of 
roughness levels. Other software will quantify the roughness characteristics of the resulting B-spline 
representations, including the characteristic distance of the curve from the data. 

Summary 

A FORTRAN 77 program that generates a smoothed B-spline representation of a given ice geometry has 
been developed. The user specifies a maximum tolerance or a root-mean-square tolerance along with 
other necessary parameters. The program returns a B-spline curve that satisfies the given tolerance. The 
software permits the rapid generation of several B-splines that satisfy a wide range of tolerance 
requirements. This approach represents a significant improvement over the current technique of 
smoothing in Smagglce 2D, in which selected control points are deleted. This program was developed for 
possible incorporation into the next version of Smagglce 2D (vl.8). 
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Any B-spline curve that approximates the ice data to within a given tolerance is not unique. Several 
additional parameters are available to the user to produce a different B-spline representation of the ice 
data. Sometimes use of these additional parameters is necessary to obtain a useful representation of the 
ice geometry. For example, for each of three different ice geometries investigated, it was found that if the 
endpoint derivatives were free and the maximum tolerance sufficiently small, the resulting approximating 
B-spline curve possessed large wiggles near an endpoint. In each of these cases, also requiring 
satisfaction of suitable endpoint derivative constraints led to an appropriate representation. 

This Year 1 report discusses the theory behind the program and illustrates its use with respect to a given 
ice geometry. The program itself appears in Appendix C. 
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Appendix A — List of Symbols 


Symbols 

d 

max 

d 

rms 

mdata 

mknot 

n 

P 

u 

u, 

u k 

C (m) 

C \u) 

N 

P 

P, 

Q k 

R 

R* 

u 

a 

^*max 

^rms 


Definition 

Maximum separation distance between B-spline and ice data; defined in Equation (16) 

Rms separation distance between B-spline and ice data; defined in Equation (17) 

Last index of ice data 
Last index of knot vector 

Last index of control points; last /-index of basis functions /V, p (u) 

Degree of B-spline curve 
Curve parameter, 0 < u < 1 
/th knot; / = 0, . . . , mknot 

kth data parameter; defined in Equation (8); k = 0 mdata 
Two-dimensional B-spline curve; functionally dependent upon curve parameter u 
Lirst derivative curve of B-spline; C \u) = dC(u)l du 

/th B-spline basis function of degree p; / = 0 n ; functionally dependent upon curve 

parameter u 

Matrix of basis functions; defined in Equation (12) or (18) 

Matrix of control points; defined in Equation (13) or (19) 

/th two-dimensional control point; i = n 
kth two-dimensional data point; k = 0, . . . , mdata 
Matrix defined in Equation (14) or (20) 

Vector defined in Equation (15) or (21) 

Knot vector; form given in Equation (2) 

Factor used to generate new maximum tolerance from relation £ max = a * <r/ max during iteration 
to satisfy rms tolerance 

Maximum tolerance 

Rms tolerance 
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Appendix B — Program Description 


Main Program 

This software is intended to be incorporated into Smagglce 2D, which has its own GUI interface. A 
simple text-based interface with users was therefore adopted for development purposes. All input and 
output is controlled by the main program. The user sequentially enters the following data: 

1 . Run number (two digits) 

2. Name of data file 

3. Are endpoint derivatives specified? If yes, either accept default values or provide values. 

4. Degree p of B-spline (default value is 3; must be 3, 4, or 5.) 

5. Initial value of n (default value is the minimum value, which is p or p + I depending on whether 
endpoint derivatives are free or fixed) 

6. Maximum tolerance (default value is £ max = l.Ox 10“ 3 ) 

7. Enter rms tolerance if desired. If so, factor a (default value: a = 0.9 ) must also be entered. 

This concludes the user’s input. 

The data file is presumed to be an ASCII file structured as follows: 

Line 1: Number of data sets (1 for just ice data or 2 

for both clean and iced airfoil data) 

Line 2: Number of points in first data set 

Lines 3 to end x-y data of first data set 

of data set 1 : 

Next line: Number of point in second data set (if 

present) 

Next line to end: x-y data of second data set 

The first data set may represent just the ice geometry or the iced airfoil. The second data set, if present, 
represents the clean airfoil. The program reads in the data file, and if necessary separates the ice data from 
the airfoil data. Data is written to appropriate file(s), and ice data is retained in memory to calculate the 
approximating B-spline curve using the iterative procedure described previously. 

The most important file exported by the main program is nurbs??.dat. It contains the data that defines the 
approximating B-spline curve. Here, 7? denotes a two-digit run number. The file has the following 
format: 

Line 1 : p (degree) 

Line 2: n (final value) 

Lines 3: mknot (final value) 

Next mknot + 1 lines: Knot vector 

Next n + 1 lines: Control points (x-y format) 
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The table below summarizes all the files exported by the main program. Following the table are brief 
statements of puipose for each subroutine included in the program file. 


File Name 

Description 

clean.dat 

Clean airfoil data (x-y format) 

cukb77.dat 

C (u k ),k = ()...., mdata {x-y format) 

ice.dat 

Ice geometry {x-y format) 

icecv77.dat 

B-spline curve {x-y format) 

log.txt 

Log file of run 

nurbs77.dat 

NURBS data file 

sum77.txt 

Summary of run 


FindSpan 

The semi-closed interval represents the ith knot span. FindSpan is a function that, given the 

curve parameter u, returns the knot span in which u is located. The value u = 1 is an exception to the 
above definition; it is assigned knot span n. 

FindSpanA 

If u represents a ID array of curve parameter values, subroutine FindSpanA returns a ID array of the 
corresponding knot span indices. 

NBasis 

Given the scalar u, its knot span index i, degree p, and knot vector U, subroutine NBasis computes the set 
of nonzero basis functions TV,. (; u ), ...,N ip (u) and returns their values in the ID array N. 

NBasisA 

Given the ID array u, a corresponding set of knot span indices, degree p, and knot vector U, subroutine 
NBasisA computes the full set of nonzero basis functions, returning them in the 2D array NA. The jth 
column of NA corresponds to the / 1 h element of u. 

Cparam 

Given the point data Q 0 , . . . , Q mdata for the ice, subroutine Cparam calculates the set of data parameters 
tf 0 , . . . , u mdata according to Equation (8), returning their values in the ID array ukb. 

Knotvec 

Subroutine Knotvec generates a knot vector U based upon Equation (10). 

NCurve 

Subroutine NCurve computes the point (x, y) = C(u ) on the B-spline curve for the scalar curve parameter 
u. 
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NCurveA 

Subroutine NCurveA computes the set of (x,y) points on the B-spline curve that correspond to a ID array 
of curve parameters. 

Rkarray 

Subroutine Rkarray computes R p . . . , R mrfflM _j corresponding to either Equation (15) or (21), depending 
upon whether the endpoint derivatives are free or fixed, respectively. 

RightHandSide 

Subroutine RightHandSide computes R according to either Equation (14) or (20), depending upon 
whether the endpoint derivatives are free or fixed, respectively. 

ABMatrix 

Subroutine ABMatrix computes the matrix N 7 N and stores it in an upper band form that is suitable for 
the LAPACK library routine DPBSV [10]. (Subroutine DPBSV solves the linear system. Equation (11), 
based upon the Cholesky method.) The matrix N is given either by Equation (12) or (18), depending upon 
whether the endpoint derivatives are free or fixed, respectively. 

SpanTest 

Subroutine SpanTest returns in the ID array nonspan the indices of all nonconforming knot spans; i.e., 
those spans that do not satisfy the tolerance requirement |C (u k ) - Q k | < £ max . 

Deviations 

Subroutine Deviations returns the maximum separation distance <7 max , the root-mean-square deviation 
d nm , and a 2D array containing the (x,y) values corresponding to C(i7 0 ),..., C(« m<tofl ). The quantities 
r/ max ar *d d nns are respectively defined by Equations (16) and (17), and in the program correspond to 
FORTRAN variables epmax and rms. 

NewU 

Subroutine NewU generates a new knot vector by inserting a new knot at the midpoint of every 
nonconforming knot span. 

Verify 

Subroutine Verify checks each knot span to insure that each contains at least one value of the set of data 
parameters u 0 ,...,u indata . 
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Appendix C — FORTRAN 77 Program 

i k 

program globalapp 

' k 

* Program Summary 

■k 

* Given a set of Q=(x,y) data points of length mdata + 1, two 

* integers n and p (such that 3 <= p <= n << mdata) , and a maximum 

* distance criterion ep (0 < ep <= 0.1, say), the program globalapp 

* employs an iterative least-squares procedure to calculate a NURBS 

* curve of degree p that approximates the data to within the nominal 

* distance ep. Optionally, the user may also specify a root-mean 

* square distance criterion (eprms) . Endpoint derivatives may be free 

* (default) or fixed. All NURBS curves generated by this program have 

* endpoints that exactly coincide with the prescribed ice data. The 

* first and last control points are equal to the respective 

* endpoints. 


Data is presumed read in a counter-clockwise direction about the 
airfoil. This is only important in referencing upper and lower 
endpoints of the ice region. 


After entry of the requisite data, the program calculates a knot 
vector, and a do while loop is entered. A matrix and right-hand 
side vector (2 columns) is then generated with control points as 
unknowns. The matrix is banded, symmetric, and positive 
definite. A special LAPACK solver that uses the efficient Cholesky 
method is then used to obtain the unknown control points. The 
resulting NURBS curve is used to determine whether the distance 
criterion eptemp (initially, the same as ep) is satisfied. If not, 
n is appropriately increased by the insertion of knots at the 
midpoints of non-converging spans, thereby creating a new knot 
vector. The do while loop is then entered again from the top, and 
the cycle is repeated until the NURBS curve satisfies the 
tolerance requirement, or the resulting matrix is singular. If at 
any stage the newly created knot vector has a span that does not 
contain a value of the parameterized curve, the whole knot vector 
is recalculated by redistribution based upon the original 
algorithm. If the redistributed knot vector also contains one or 
more spans without a value of the parameterized curve, the do 
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while loop is exited, and the approximation fails. 


If eprms is specified, the iteration procedure is slightly 
modified as follows. The specified maximum distance specification 
is first met as indicated above. If eprms is specified, the 
root mean squared distance is compared with eprms. If this 
specification is satisfied, the iteration loop is exited. If not, 
the maximum distance specification eptem is reduced by a factor alph 
0 < alpha < 1, and iterations continue until this new maximum 
distance specification is met. The rms specification is then 
checked again. The cylce repeats until the rms specification 
is satisfied. 


~ k 
k 
k 
k 
k 
k 
k 
k 
k 
k 
k 
k 
k 
k 
k 
k 
k 
k 
k 
k 
k 
k 
k 


The procedure when endpoint derivatives are specified is nearly 
the same. An endpoint derivative is fully specified by a polar 
angle and a magnitude. By default, first-order finite difference 
calculations are used to prescribe these derivatives. After 
a knot vector is first defined as above, the derivative 
information is used to determine the 2nd and next to last 
control points. The do while loop is then entered as above. 

Though governing equations differ slightly in this case, a 
linear system of equations is solved to find the remaining 
unknown control points. A modified knot vector is then created 
in a manner similar to that above. If at any time the knot 
vector needs to be redistributed, the second and second to last 
control points are recalculated. The fixed derivative algorithm 
was specially developed for this program, and is not found in the 
reference below. 

After a NURBS curve is found that satisfies specified tolerance ( s ) , 
the requisite information of the NURBS curve is saved 
to a file for later reconstitution. A log file is also kept, 
as well as a file that summarizes the results of the analysis. 

Also output is an x-y data file that represents the generated 
NURBS curve for convenient graphing. 




* Reference 

k 

* Piegl, L. and Tiller, W., The NURBS Book, 2nd edition. New York: 

* Springer-Verlag, 1997. 

k 
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* NOTE: In Piegl and Tiller, the first index of most vectors begin 

* with zero. To avoid confusion, we follow that convention as 

* appropriate. We also use indices 1 and 2 to denote x and y 

* components of certain arrays. 

* 

ici(i(i(i(i(ici(ici('k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 


~k 
■k 
: k 
' k 
■k 
■k 
■k 
: k 
k 
k 


Record of Revisions: 


Date 


Programer 


05/20/03 Loren H. Dill 

05/22/03 L.H. Dill 

08/06/03 L.H. Dill 


Description of Change 


Original code 
Added rms tolerance 
Increased saved digits 
in output (nurbs??.dat) 


k k 


* Parameters 


~k 


i k 
' k 
■k 
: k 
' k 
' k 
■k 
■k 
■k 
: k 


mcleandatamax 

mdatamax 

m2datamax 

nmax 

nredistmax 

pmax 


Largest index of x-y points output to data file 
for clean airfoil (file: clean.dat) 

Largest anticipated index of input data vector, 
clean plus ice 

Largest index of output x-y data vector for rough 
ice (along NURBS curve) (file: icecv??.dat) 

Largest anticipated index of control points 
Largest index of nvalues vector 
Highest anticipated degree of NURBS curve 


* Major Input Variables 

i k 


mdata Number of x-y data points (less 1 after initial index set 
to 0) . mdata is initially read in, and depending upon 
dataset format, may be total number of ice+clean data 
points. Thoughout the main part of the program, mdata 
represents the number of ice data less one. 

Q x,y ice data points — a 2D array 

p Desired degree of NURBS curve 

n Largest index of NURBS curve. Initially given and then 

increased during the iteration procedure 
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■k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 


phiO, phiM First and last polar angles corresponding to desired 

slopes of NURBS curve at beginning and end. (in degrees) . 

DOL, DmL Magnitude of first derivative vectors of NURBS curve 
at beginning and end 

sDOL, sDmL Scale factors, relative to default values, of the 

magnitude of first derivative vectors of NURBS curve 
at beginning and end 

ep Desired maximum tolerance between NURBS curve and prescribed 
data 

eprms Desired maximum rms tolerance between NURBS curve and 
prescribed data 

run A two-digit character array designated the run number. Used 
to associate output file names with given runs. 

infilename A character array of length 12 used to identify 
the file containing the x-y ice data. 




k 


* Major Working Variables 

k 


k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 


epmax Calculated value of maximum distance from curve to the 
ice data. Actually a nominal value, 
eptemp If only maximum tolerance is specified, eptemp = ep . 

If eprms tolerance is specified, eptemp < ep if eprms 
is not achieved. eptemp is periodically reduced until 
eprms is achieved. 

mknot Largest index of knot vector U, mknot = n+p+1 

INFO Variable from LAPACK routine dpbsv that solution 

of linear matrix equation was successful or not 
inonspan number of nonconverging knot spans 
nonspan vector containing spans that have not converged 
nvalues a 1-D array containing the n values for which 
the knot vector was redistributed. Only last 
nredistmax values are retained, 
offset Derivative flag. 

offset = 1 if endpoint derivatives are not set. 
offset = 2 if they are 
P Control points 

ukb array containing values of parameterized curve; each 

value of array corresponds to respective data point 
spanA an array that contains the span indices of ukb relative 
to a given knot vector, 
u Curve parameter 

U Knot vector 
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■k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 


Unew A new knot vector awaiting verification before acceptance 
NA Array containing the non-zero basis functions for each 
element of a parameter array 

Rk array containing the Rk values of Piegl & Tiller, p 411, 
if end derivatives are free, and a modified version 
if end derivatives are fixed. 

R array representing right-hand side of eq. 9.65 of Piegl 

& Tiller, p 411, if end derivatives are free, and a 
modified version if end derivatives are fixed. 

NTNB matrix (in banded form) for linear eq. 9.65 of P&T, 
if end derivatives are free, and a modified version 
if end derivatives are fixed. 

CA array containing the x-y data of the NURBS curve 
corresponding to a given parameter array 
DO, DM Vectors representing endpoint derivatives if prescribed 
rms the root mean square distance between the curve and the 

set of data points. Again a nominal value, 
success A logical variable indicating that a newly created knot 
vector has an element of ukb within each span, or not. 




k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 

k 


OUTPUT FILE DESCRIPTIONS 


clean . dat 
cukb? ? . dat 

ice . dat 
icecv? ? . dat 
log . txt 
nurbs?? .dat 
sum? ? . txt 


The clean airfoil data, if given upon input 
Data file containing x-y data of NURBS curve 
corresponding to paramter array ukb 
The ice data 

The NURBS curve x-y data 

A log file for the last execution of glapp 
A file containing NURBS data, such as the knot vector 
A file that summarizes a run 


In the above, ?? denotes the run number given upon input. 


kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

* Compile Command for g77 compiler 

*g77 -f source-case-preserve -Wunused -Wall -Wsurprising glapp. f -llapack 
lblas -o glapp 

* Use xmgrace or other plotting program to examine output *dat files 


implicit none 


integer mdata, mdatamax, m2datamax, nmax, pmax, mcleandatamax 
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integer ndatasets, nredist, nredistmax, mcleandata, nmaxm 

parameter (mdatamax=600 , nmax=600 , pmax=5 , m2datamax=8 1 92 ) 
parameter (mcleandatamax=250 , nredistmax=5 ) 

integer i , j , mknot , of f set , p, n, spanA ( 0 :m2datamax) 
integer nonspan (nmax-1 ), inonspan 
integer nvalues (nredistmax) , nf irst, nlast 
integer INFO 

double precision U ( 0 : nmax+pmax+1 ) , du, u ( 0 :m2datamax) 

double precision Q (2, 0 rmdatamax) , Qclean (2, 0 :mcleandatamax) 

double precision ukb ( 0 rmdatamax) , P (2 , 0 : nmax) , Unew ( 0 : nmax+pmax+1 ) 

double precision NA (0 rpmax, 0 rmdatamax) 

double precision CA (2, 0 rm2datamax) 

double precision Rk (2, 0 rmdatamax) , R (nmax-1, 2) 

double precision NTNB (pmax+1 , nmax-1 ), ep, eprms, eptemp, alpha 

double precision DO (2 ) , Dm (2 ) , sDOL, sDmL 

double precision phiO, phim, DOL, DmL, pi , epmax, rms 

character*12 infilename 

character*l ans, ansD, ansDC, ansDr, ansphi, ansdeg, ansn, ansrms 
character*2 run 

logical success 

Open all output files and write a blank to clear data from 
any previous run of same number. Open and close later as needed, 
write (*,*)' Please enter a 2-digit run number.' 
read ( * , * ) run 

open ( 7 , f ile= ' icecv ' //run/ / ' . dat ' ) 

open ( 8 , f ile= ' log . txt ' ) 

open ( 9, f ile= ' clean . dat ' ) 

open (11, f ile= ' ice . dat ' ) 

open ( 12 , f ile= ' nurbs ' / / run/ / ' . dat ' ) 

open (13, f ile= ' sum ' / /run/ / ' . txt ' ) 

open ( 14 , f ile= ' cukb ' / /run/ / ' . dat ' ) 

write ( 7 , * ) ' ' ; write ( 8 , * ) ' ' ; write ( 9, * ) ' ' 

write (11,*) ' ' ; write ( 12 , * ) ' ' ; write ( 13, * ) ' ' 

write ( 14 , * ) ' ' 

close (7) ; close (8) ; close (9) 

close (11) ; close (12) ; close (13) ; close (14) 
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Read in airfoil data. First line of data file contains number of 
datasets (1 = only ice data; 2 = both clean and ice data) . 

Second line contains number of data points of first data set. Ask 
user for name of input file then read first two lines. 

write (*,*)' Enter name of data file as string.' 
read ( * , * ) infilename 

open (10, f ile=inf ilename) 
read (10, *) ndatasets 
read ( 10 , * ) mdata 

mdata = mdata -1 ! Adjust maximum because indices begin at 0 

if (ndatasets .eq. 1) then 

if (mdata . gt. mdatamax) then 

write (*,*) 'Number of data points ', mdata, ' exceeds maximum 

1 , 'expected ', mdatamax, ' . Aborting . . . ' 

open ( 13 , f ile= ' sum ' / /run// ' . txt ' ) 

write (13,*) 'Number of data points ', mdata, ' exceeds maximum 

1 , 'expected ', mdatamax, ' . Aborting . . . ' 

close (13) 
stop 
end if 

do i=0, mdata 

read (10, *)Q(l,i),Q(2,i) 
end do 
close ( 10 ) 

Write ice data to file 

open ( 1 1 , f ile= ' ice . dat ' ) 
do i=0, mdata 

write (11, 10) Q (1, i) ,Q(2,i) 
end do 
close (11) 

else if (ndatasets .eq. 2) then 
mcleandata = mdata 

if (mcleandata .gt. mcleandatamax -1) then 

write (*,*) 'Number of clean airfoil data points ', mcleandata 

1 , ' exceeds maximim expected ', mcleandatamax, ' . 

2 ' Aborting . . . ' 

open ( 13, f ile= ' sum ' / /run/ / ' . txt ' ) 

write (13,*) 'Number of clean airfoil data points ',mcleandat 

1 , ' exceeds maximim expected ', mcleandatamax, ' . ', 

2 ' Aborting . . . ' 
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close (13) 
stop 
end if 

do i=0 , racleandata 

read (10, *) Qclean (1, i) , Qclean (2, i) 
end do 

open (9, file= 'clean. dat ' ) 
do i=0 , racleandata 

write (9, 10) Qclean (1, i) , Qclean (2, i) 
end do 
close ( 9 ) 
read ( 10 , * ) mdata 

if (mdata .gt. mdatamax-1) then 

write (*,*) 'Number of data points ', mdata, ' exceeds maximum ' 
1 , 'Increase mdatamax. Aborting . . . ' 

open ( 13, f ile= ' sum ' / /run/ / ' . txt ' ) 

write (13,*) 'Number of data points ', mdata, ' exceeds maximum' 
1 , ' Increase mdatamax. Aborting . . . ' 

close (13) 
stop 
end if 

mdata = mdata - 1 
do i=0, mdata 

read (10, *)Q(l,i),Q(2,i) 
end do 
close ( 10 ) 

Check that the length of the ice data exceeds the length 
of the clean airfoil data by a sufficient amount. The 20 
is not a rigid requirement. This is only to separate ice and 
clean airfoil data 

if (mdata .it. mcleandata+20 ) then 

write (*,*)' Ice data length too short to proceed.' 
write (*,*) 'Aborting' 
stop 
end if 

Determine where ice begins in ice/airfoil data 
i=0 

do while ( (abs (Qclean (1, i) -Q (1, i) ) .it. 1. e-6 ) .and. 

1 (abs (Qclean (2, i) -Q (2, i) ) .it. 1. e-6 ) ) 

i=i + l 
end do 
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nf irst=i 


Determine where ice ends in ice/airfoil data 
i=0 

do while ( (abs (Qclean ( 1 , mcleandata-i ) -Q (l,mdata-i) ) . it . 1 .e-6 ) 

1 . and . (abs (Qclean (2 , mcleandata-i ) -Q (2 , mdata-i )). it . 1 . e-6 )) 

i=i + l 
end do 
nlast=i 

mdata = mdata - nfirst - nlast 

Transfer ice data to beginning of arrays 
do i = 0, mdata 

Q(l,i)=Q(l,nfirst+i) 

Q(2,i)=Q(2,nfirst+i) 
end do 


Write ice data to file 

open ( 1 1 , f ile= ' ice . dat ' ) 
do i=0, mdata 

write (ll,10)Q(l,i) ,Q(2,i) 
end do 
close (11) 
end if 

Begin summary and log files. Get input file name and run number. 

open (13, f ile= ' sum ' / /run/ / ' . txt ' ) 

write (13,*) 'Input data filename = ' , inf ilename 

write (13,*) 'Run number = ' , run 

open ( 8 , f ile= ' log . txt ' ) 

Parameterize data 

call Cparam (mdata, Q, ukb) 


Inquire about endpoint derivatives 
pi=4 . 0 *atan (1.0) 
ans = 'O' 

do while (ans .ne. '1') 

write (*,*) 'Do you want endpoint derivatives specified (y/n)?' 
read ( * , * ) ansD 

if (ansD .eq. 'Y' .or. ansD .eq. 'y')then 

DO ( 1 ) = (Q ( 1 , 1 ) -Q ( 1 , 0 ) ) /ukb ( 1 ) ! Calculate default values 

DO (2) = (Q (2, 1) -Q (2, 0) ) /ukb (1) 

Dm ( 1 ) = (Q ( 1 , mdata) -Q ( 1 , mdata-1 ) ) / ( 1-ukb (mdata-1 ) ) 
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Dm (2 ) = (Q (2 , mdata) -Q (2 , mdata-1 ) ) / ( 1-ukb (mdata-1 ) ) 

DOL = sqrt (DO (1) *D0 (1) +D0 (2) *D0 (2) ) 

DmL = sqrt (Dm ( 1 ) *Dm ( 1 ) +Dm (2 ) *Dm (2 ) ) 
phi0=atan2 (DO (2) , DO (1) ) *180. /pi 
phim=atan2 (Dm (2 ) , Dm (1) ) *180. /pi 
write (13, * ) 

write (13,*) 'Endpoint derivatives specified.' 

write (*,*) 'Do you want endpoint derivatives calculated ' 

write (*,*) 'by first-order differencing of first and last ' 

write (*,*) 'endpoint data pairs (y/n)?' 

read ( * , * ) ansDC 

if (ansDC .eq. 'Y' .or. ansDC .eq. 'y')then 

write (13,*) 'Derivatives calculated by differencing' 
write (13,*) 'of endpoint data pairs' 
ans = ' 1 ' 
of f set=2 

else if (ansDC .eq. 'N' .or. ansDC .eq. 'n') then 
write ( * , * ) 

write (*,*) 'Both direction and magnitudes of endpoint' 
write (*,*) 'derivative vectors must be specified. Do' 
write (*,*) 'you want the default directions determined' 
write (*,*) 'from finite differencing of the first and' 
write (*,*) 'last pairs of prescibed data? (y/n) ' 
read ( * , * ) ansphi 

if (ansphi .eq. 'N' .or. ansphi .eq. 'n') then 
write (*,*) 'Default values for upper and lower' 
write (*,*)' angles in degrees:' 
write ( * , 50 ) phiO , phim 
write ( * , * ) 

write (*,*) 'Enter upper and lower endpoint' 
write (*,*)' directions as polar angles in degrees:' 
read ( * , * ) phiO , phim 
end if 
write ( * , * ) 

write (*,*) 'Do you want finite differences of endpoint' 
write (*,*) 'pairs to determine the magnitudes of ' 
write (*,*) 'these derivatives (y/n)?' 
read ( * , * ) ansDr 

if (ansDr .eq. 'N' .or. ansDr .eq. 'n') then 
write ( * , * ) 

write (*,*) 'Specify the upper and lower ' 
write (*,*) 'magnitudes of these vectors as' 
write (*,*)' scale factors of the default values:' 
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read ( * , * ) sDOL, sDmL 
DOL=sDOL*DOL 
DmL= s DmL * DmL 
end if 

DO (l)=D0L*cos (phi0*pi/180 . ) 

DO (2)=D0L*sin (phi0*pi/180 . ) 

Dm (1) =DmL*cos (phim*pi/180 . ) 

Dm (2 ) =DmL*sin (phim*pi/ 180.) 

write (13,*) 'Derivatives specified by user: 

ans = ' 1 ' 

of f set=2 


end if 



write (13, * ) 

'Polar angle (deg 

write (13, * ) 

' upper: 

' , phiO 

write (13, * ) 

' lower: 

' , phim 

write (13, * ) 

' Lengths : 

f 

write (13, * ) 

' upper: 

' , DOL 

write (13, * ) 

' lower: 

' , DmL 

write (13, * ) 




else if (ansD .eq. 'N' .or. ansD .eq. 'n') then 
write (13, *) 

write(13,*) 'Endpoint derivatives not specified by user.' 
write (13, *) 
ans = ' 1 ' 
offset = 1 
end if 
end do 

Set degree p and initial number of terms via n 
p=3 

write (*,*) 'The default degree of the NURBS curve is 3.' 
write (*,*) 'Do you want the default degree (y/n) ' 
read ( * , * ) ansdeg 

if (ansdeg .eq. 'N' .or. ansdeg .eq. 'n') then 

write (*,*) 'Degree must satisfy 3 <= p <= pmax = ' , pmax 
write ( * , * ) ' Enter degree : ' 
read ( * , * ) p 

if ( p . gt . pmax .or. p .It. 3) then 

write (*,*) 'NURBS degree p = ',p, ' not in range' 
write ( * , * ) ' Aborting run . . . ' 

write (13,*) 'NURBS degree p = ',p, ' not in range' 
write ( 13 , * ) ' Aborting run . . . ' 
close (13) 
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stop 
end if 
end if 

nmaxm = min (nmax, mdata - p - 2) 

if ( (offset .eq. 1 .and. nmaxm .It. p) .or. 

1 (offset .eq. 2 .and. nmaxm .It. p+1) ) then 

write (*,*)' Your dataset is too small to proceed. Aborting..' 
write ( 8 ,*)' Your dataset is too small to proceed. Aborting..' 
write ( 13, *)' Your dataset is too small to proceed. Aborting..' 
stop 
end if 
write ( * , * ) 

write (*,*) 'For the number of terms in the initial NURBS' 

write (*,*) 'curve, do you want the minimum acceptable value (y/n)?' 

read ( * , * ) ansn 

if (ansn .eq. 'Y' .or. ansn .eq. 'y')then 
if (offset .eq. l)n=p 
if (offset .eq. 2)n=p+l 
else 

write (*,*) 'Enter n to specify number of terms in ' 

1 , 'initial NURBS curve' ! n actually specifies number 

if (offset .eq.l) then ! terms less one 

write (*,*) 'Initial n must satisfy ',p, ' < = n < ', nmaxm 

else 

write (*,*) 'Initial n must satisfy ',p+l, ' < = n < ', nmaxm 

end if 
write ( * , * ) 

write (*,*) ' (This maximum is an upper limit. Typically want' 
write (*,*)' to set n near lower limit and n << ', nmaxm, '. ' 
write (*,*) 'Matrix may be singular for n < ', nmaxm, ' depending ' 
1 , 'upon data set.) ' 

write ( * , * ) 

write (*,*) 'What value of n do you want?' 
read ( * , * ) n 

if ( (offset .eq. 1 .and. n .It. p) .or. 

1 (offset .eq. 2 .and. n .It. p+1) .or. 

2 n . gt . nmaxm ) then 

write (*,*)'n is outside specified range.' 
write ( * , * ) ' Aborting run . . . ' 

write ( 13, * ) ' Specified n = ' , n, ' is outside range . ' 
write (13, *) 'Aborting run . . . ' 
close (13) 
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stop 
end if 
end if 

write(13,*) 'Initial n: ',n 

Inquire about desired maximum tolerance 
ans = 'O' 

do while (ans . ne. '1') 

write (*,*) 'Do you want the default maximum tolerance ', 

1 ' (ep = 0.001) (y/n) ?' 

read ( * , * ) ans 

if (ans .eq. 'Y' .or. ans .eq. 'y')then 
ep= 1.0d-3 
ans = ' 1 ' 

else if (ans .eq. 'N' .or. ans .eq. 'n') then 

write (*,*) 'Enter maximum tolerance (0 < ep <= 0.1) ' 
read ( * , * ) ep 
ans = ' 1 ' 
end if 
end do 
write ( * , * ) 

if ( (ep .gt. 0.1 ) .or. (ep .le. 0.) ) then 

write (*,*) 'Distance criterion ep = ',ep, ' should be ' 

1 , ' in range 0 < ep <= 0.1. Aborting . . . ' 

write (13,*) 'Distance criterion ep = ',ep, ' should be ' 

1 , ' in range 0 < ep <= 0.1. Aborting . . . ' 

close (13) 
stop 
end if 

Set eptemp equal to ep . eptemp will only differ from ep 
if root-mean square tolerance is set. 
eptemp = ep 

Inquire about root-mean square tolerance 
ans = 'O' 

do while (ans .ne. '1') 

write (*,*) 'Do you want to specify a root mean square tolerance ', 
1 ' (y/n) ? ' 

read ( * , * ) ansrms 

if (ansrms .eq. 'Y' .or. ansrms .eq. 'y')then 

write (*,*) 'Enter the root mean square tolerance. Your value ', 
1 'should lie between 0 and ' , ep 

read ( * , * ) eprms 

if (eprms . gt . 0. .and. eprms .It. ep)then 
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alpha = 0.9 

write (*,*) 'The default factor by which to reduce the' 
write (*,*) 'maximum deviation when attempting to satisfy' 
write (*,*)' your root mean square tolerance is 0.9.' 
write (*,*) 'Do you accept the default factor?', 

1 ' (y/n)?' 

read ( * , * ) ans 

if (ans .eq. 'N' .or. ans .eq. ' n ' ) then 
write (*,*)' Enter the desired factor.' 
read ( * , * ) alpha 

if (alpha .gt. 0. .and. alpha .It. l)then 
ans = ' 1 ' 
end if 
else 

ans = ' 1 ' 
end if 
end if 

write ( 8 ,*)' Reduction factor to achieve eprms alpha 
else 

ansrms = 'n' 
ans = ' 1 ' 
end if 
end do 

write(8,*) 'p , n, ep = ',p, ' ',n, ' ' , ep 

if (ansrms .eq. 'Y' .or. ansrms .eq. 'y')then 
write (8,*) 'eprms = ', eprms 
end if 

write(8,*) 'x,y, ukb data' 
do i = 0,mdata 

write(8,10)Q(l,i) ,Q(2,i) , ukb (i) 
end do 

write (13,*) 'Specified NURBS degree: ', p 

write (13,*) 'Specified nominal maximum tolerance: ' , ep 

if (ansrms .eq. 'Y' .or. ansrms .eq. 'y')then 
write (13,*) 'Specified rms tolerance: ', eprms 

end if 


Create knotvector based upon Piegl & Tiller, p. 412, eqn. (9.69) 
and calculate mknot (number of knots less one) . Keep track of 
value of n when Knotvec is called via vector nvalues. 
nredist=l 
nvalues ( 1 ) =n 
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call Knotvec (mdata,n, p, ukb, U) 
mknot = n + p + 1 
write ( 8 , * ) 

write (8,*) 'Knot Vector U(i) ' 
write(8,20) ( i , U ( i ), i=0 , mknot ) 
write ( 8 , * ) 

* Verify each knot span contains at least one value of ukb. Output 

* argument success returns .true, if all is OK. If unsuccessful, 

* execution is stopped. 

call Verify (mdata, n,p,ukb,U, success) 

if (success . eqv. .false.) then 

write (8,*) 'Before entering main loop, not all knot spans', 

1 'contain at least one value of ukb. Aborting . . . ' 

write (*,*) 'Before entering main loop, not all knot spans', 

1 'contain at least one value of ukb. Need to reduce', 

2 'initial value of n. Aborting . . . ' 

write (13,*) 'Before entering main loop, not all knot spans', 

1 'contain at least one value of ukb. Need to reduce', 

2 'initial value of n. Aborting . . . ' 

close (13) 

stop 
end if 

* Enter iteration loop to calculate global approximating NURBS 

* curve. Require at least one iteration by setting number of non- 

* converging spans greater than zero, e.g., inonspan = 1. Loop 

* repeats until the number of non-converging spans drops to zero 

* (inonspan =0), at least one knot span in new knot vector does not 

* contain a parameter value, or the number of terms in the next 

* NURBS representation exceeds nmax + 1 (i.e., n > nmax ) . 

inonspan = 1 

write (8,*) 'Entering do while' ! Main iteration loop, 
do while ( inonspan . gt . 0 .and. success .and. n .le. nmax ) 
write ( 8 , * ) 

write (8,*) 'Top of do while ' 

write ( 8 , * ) ' n = ',n,' mknot = ' , mknot 
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First, find span indices of all elements of ukb vector. 

Next, calculate all non-zero basis functions associated with 
elements of ukb. 


call FindSpanA (mdata, n, p, ukb, U, spanA) 

call NBasisA (mdata, mknot , p, pmax, spanA, ukb, U, NA) 


■k 
: k 
k 


k 


The first and last control points are simply the specified 
endpoint data. If end derivatives are specified, the second 
and next to last control points are determined by the 
derivative information. 


P (1,0) =Q (1,0) 

P (2, 0) =Q (2,0) 
if (offset .eq. 2) then 

P ( 1 , 1 ) =Q (1,0)+ U (p+1 ) / p * DO ( 1 ) 

P (2, 1) =Q (2, 0) + U (p+1 ) / p * DO (2 ) 

P (l,n-l) =Q (1, mdata) - (1.0 - U(n))/p * Dm ( 1 ) 
P (2, n-1) =Q (2, mdata) - (1.0 - U(n))/p * Dm(2) 
end if 

P (1, n) =Q (1, mdata) 

P (2, n) =Q (2, mdata) 


Next, calculate the set of Rk arrays, Piegl & Tiller, eqn. 
(9.63) on page 411, and then the right-hand side of (9.65) : 
call Rkarray (mdata, n, of f set, p, pmax, spanA, Q, ukb, P, U, NA, Rk) 
call RightHandSide (mdata, n, nmax, of f set, p, pmax, spanA, NA, Rk, R ) 

Calculate the NTNB matrix of eqn. (9.65) in banded form. 

call ABMatrix (mdata, n, offset, p, pmax, spanA, NA, NTNB) 

Call the LAPACK banded matrix solver dpbsv to find the 
locations of the internal control points. (The first and last 
control points are simply the first and last data points. 

If endpoint derivatives are specified, the next inner control 
points are determined from derivative information.) . If 
the solver fails (INFO .ne. 0), an error message is generated. 

call dpbsv (' u ', n-2 *of f set+1 , p, 2, NTNB, pmax+1, R, nmax-1, INFO ) 

if (INFO .eq. 0)then 
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do j=of f set, n-of f set 

P (1, j ) =R ( j+l-of f set , 1) 

P (2 , j ) =R ( j+l-of f set , 2 ) 
end do 

else if (INFO .It. 0)then 

write (8,*) 'The ',-INFO, ' argument had an illegal value', 

1 ' in dpvsv. Aborting . . . ' 

write (*,*) 'The ',-INFO, ' argument had an illegal value', 

1 ' in dpvsv. Aborting . . . ' 

write (13,*) 'The ',-INFO, ' argument had an illegal value', 
1 ' in dpvsv. Aborting . . . ' 

close (13) 
stop 
else 

write (8,*) 'Matrix in dpbsv is singular (U (', INFO, INFO, 

1 ' ) = 0 . Aborting . . . ' 

write (*,*) 'Matrix in dpbsv is singular (U (', INFO, INFO, 

1 ' ) = 0 . ' 

write (*,*) 'May need to increase tolerance ep or ' 
write (*,*)' decrease initial n. Aborting ...' 

write (13,*) 'Matrix in dpbsv is singular (U (', INFO, INFO, 

1 ' ) = 0 . ' 

write (13,*) 'May need to increase tolerance ep or ' 
write ( 13, *)' decrease initial n. Aborting ...' 
close (13) 
stop 
end if 
write ( 8 , * ) 

write(8,*) 'Control points P(l, j) & P(2, j) for n = ',n 
write (8, 30) (j,P(l,j),P(2,j),j=0,n) 


The distance |C( ukb(k) ) - Q(k) | for each k = 1, mdata -1 is 
compared with distance specification ep (or eptemp if attempting 
to satisfy rms requirement) . Here, Q(k) represents 
the kth data point of the airfoil. If this distance exceeds the 
specification for any data point within a given span, the entire 
span is declared to be non-converging and is tagged for 
subdivision for the next iteration. 

call SpanTest (eptemp, mdata, n,p,Q,ukb,P,U, non span. 
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If all spans have converged (inonspan = 0), we calculate maximum 
(epmax) and root mean square (rms) deviations between curve and 
data. If eprms tolerance is not specified, the do while loop is 
exited. Otherwise, calculated and specified values for the 
root mean square error is compared. If rms > eprms, eptemp 
is set equal to alpha * epmax and SpanTest is called again. 

Here, alpha is a give parameter that lies between 0 and 1. Thus, 
we are guaranteed that at least one span will be non-compliant ; 
i.e., inonspan will be greater than 0 as determined by SpanTest. 

Note: If alpha is set too low, the tolerance eptemp may become 

unnecessarily small, and the program could fail. On the other 
hand, if alpha is too close to unity, eptemp will be reduced by 
only a small amount, and most likely only one span will be 
non-compliant. Hence, lots of iterations of the do while 
loop might be necessary for convergence to eprms. The default 
value of alpha = 0.9 should insure rapid convergence and in 
most cases eptemp should not become so small as to cause a 
singular matrix. If the procedure does bomb out, either eprms 
or alpha needs to be increased appropriately. 

If inonspan > 0, we determine whether n will exceed nmax if 
another iteration is performed. If it will, execution is stopped 
If the new n will be less than nmax, a new knot vector (Unew) 
is calculated. The new knot vector is formed by adding a new knot 
to the midpoint of each non-converging span. This new vector 
is tested to insure each span contains at least one value of 
ukb. If it does, the do while loop is executed again. Otherwise 
a new knot vector is created for which all knots are 
redistributed. If this new knot vector contains any spans 
lacking a value of ukb, execution is stopped because the matrix 
NTNB is no longer guaranteed to be positive definite and well 
conditioned . 

if (inonspan .eq. 0) then 
if (eptemp .eq. ep)then 
write ( 8 , * ) 

write (8,*) 'Maximum tolerance achieved.' 
write (13, * ) 

write (13,*) 'Maximum tolerance achieved.' 
end if 
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call Deviations (mdata, n, p, pmax, Q, ukb, P, U, epmax, rms, CA) 
if ( (ansrms .eq. 'Y' .or. ansrms .eq. 'y') .and. 

1 rms . gt . eprms) then 

write ( 8 , * ) 

write ( 8 ,*)' epmax, rms = epmax, ' ' , rms 

eptemp = alpha * epmax 
write (8,*) 'eptemp = ', eptemp 

call SpanTest (eptemp, mdata, n,p,Q,ukb,P,U, non span, 

1 inonspan) 

else if ( (ansrms .eq. ' Y' .or. ansrms .eq. 'y') .and. 

1 rms .le. eprms) then 

write (8,*) 'rms tolerance achieved' 
write ( 8 ,*)' epmax, rms = ', epmax, ' ' , rms 

end if 
end if 

if there are non-converging spans, and new n > nmax: 
if ( inonspan . gt . 0 .and. n + inonspan . gt . nmax ) then 

write (8,*) 'Global approximation not converged, and next ' 

1 , 'iteration will exceed maximum n. Aborting. . . ' 

write (13,*) 'Global approximation not converged, and next ' 

1 , 'iteration will exceed maximum n. Aborting. . . ' 

close(13) ;close(8) 
stop 

if there are non-converging spans, and new n <= nmax: 
elseif ( inonspan . gt . 0 .and. n + inonspan .le. nmax ) then 
call NewU (mdata, n, non span, inonspan, p, ukb, U, Unew) 
mknot = n+p+1 
write ( 8 , * ) 

write (8,*) 'Potential new U for n = ',n 
write (8,*) 'New mknot = ', mknot 
write (8, 20) (i, Unew (i) , i=0, mknot) 
write ( 8 , * ) 

call Verify (mdata, n, p, ukb, Unew, success) 

if (success) then ! The new knot vector created by inserting 

do i=0, mknot ! new knots at midpoints of non-converging 

U(i)=Unew(i) ! spans is OK. Save new knot vector 
end do ! to log file, 

write ( 8 , * ) 

write (8,*) 'New Knot Vector OK for n = ',n 
else 

write (8,*) 'Modified knot vector had span with no ', 

1 'parameter value. Redistributing knots to ', 
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2 


'create all new knot vector. ' 


i k 


~k 


nredist=nredist+l 

if (nredist . le. nredistmax) then 
nvalues ( nredist ) = n 
else 

do i = 1, nredistmax -1 

nvalues ( i ) =nvalues ( i+1 ) 
end do 

nvalues (nredistmax) = n 
end if 

call Knotvec (mdata,n, p, ukb, Unew) 
write ( 8 , * ) 

write(8,*) 'All new knot vector Unew(i) for n = ',n 

write(8,*) 'Before Verify' 

write (8, 20) (i, Unew (i) , i=0,mknot) 

write ( 8 , * ) 

call Verify (mdata, n, p, ukb, Unew, success) 
if (success) then ! Redistributed knot vector is OK. 
do i=0,mknot 
U ( i ) =Unew ( i ) 
end do 

write (8,*) 'New Knot Vector OK for n = ',n 
If endpoint derivatives are specified, define new control 
pts P(l) and P(n-l) 

if (offset .eq. 2) then 

P ( 1 , 1 ) =Q (1,0)+ U(p+l)/p * DO ( 1 ) 

P (2, 1) =Q (2, 0) + U(p+l)/p * DO (2 ) 

P (1, n-1) =Q (1, mdata) - (1.0 - U(n))/p * Dm ( 1 ) 

P (2, n-1) =Q (2, mdata) - (1.0 - U(n))/p * Dm(2) 
end if 
else 


write(8,*) 'Totally new knot vector still has', 

1 ' span with no ukb value. ' 

write(*,*) 'Totally new knot vector still has', 

1 ' span with no ukb value. ' 

write(13,*) 'Totally new knot vector still has', 

1 ' span with no ukb value. ' 

if (eptemp .eq. ep) then 

write (8,*) 'Maximum tolerance specification ', 

1 ep, ' is too small' 

write (13,*) 'Maximum tolerance specification ', 
1 ep, ' is too small' 

else 
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1 


write (8,*) 'RMS error spec is too small. Actually', 

' achieved rms of ' , rms 

write (13,*) 'RMS error spec is too small. Actually', 
1 ' achieved rms of ' , rms 

end if 
close (13) 
stop 
end if 
end if 
end if 
end do 

write ( 8 , * ) ' n, mknot ep = ' , n, ' ' , mknot, ' ' , ep 

write ( 8 , * ) 

write(8,*) 'Control points P(l, j) & P(2, j) for n = ',n 
write (8,30) (j,P(l,j),P(2,j),j=0,n) 
write ( 8 , * ) 

write (8,*) 'Knot Vector U(i) ' 
write(8,20) ( i , U ( i ), i=0 , mknot ) 
close ( 8 ) 

open ( 12 , f ile= ' nurbs ' / / run/ / ' . dat ' ) 

write ( 12 , 30 ) p 

write ( 12 , 30 ) n 

write ( 12 , 30 ) mknot 

do i = 0, mknot 

write (12, 10) U (i) 
end do 
do i = 0 , n 

write(12,10)P(l,i) , P (2, i) 
end do 
close ( 12 ) 

if (ansrms .eq. 'Y' .or. ansrms .eq.'y') then 
write (13,*) 'RMS tolerance spec achieved.' 
end if 

write (13,*) 'Knot vector redistributed for ' 
write(13,*) 'following values of n. If redistribution' 
write (13,*) 'occured more than ' , nredistmax, ' times, ' 
write (13, *)' only last ', nredistmax, ' times are reported.' 
write (13, 40) (nvalues (i) , i=l,min (nredist, nredistmax) ) 
write (13,*) 'Knot vector was redistributed a total of ', nredist, 
1 ' times . ' 

write (13, *) 

write(13,*) 'Final value of n: ',n 
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write (13, * ) 

if (offset .eq. 1) then 


DO (1) 

=p* (P 

(1 

,D-P(1, 

0 

) ) /U(p+1) 

DO (2) 

=p* (P 

(2 

,1)-P(2, 

0 

) ) /U(p+1) 

Dm ( 1 ) 

=p* (P 

(1 

, n) -P ( 1 , 

n 

-1) ) / (1.0-U(n) ) 

Dm (2 ) 

=p* (P 

(2 

, n) -P (2, 

n 

-1) ) / (1.0-U(n) ) 

D0L = 

sqrt 

(DO (1) *D0 ( 

1 

) +D0 (2) *D0 (2) ) 

DmL = 

sqrt 

(Dm ( 1 ) *Dm ( 

1 

) +Dm (2 ) *Dm (2 ) ) 

phi0 = 

=atan2 

(DO (2) , DO 

1 

) ) *180 . /pi 

phim= 

=atan2 

(Dm (2 ) , Dm ( 

1 

) ) *180 . /pi 

write ( 13, * 

) ' 

Endpoint 

Derivative Information 

write ( 13, * 

) ' 

Polar an 

gle (degrees) : ' 

write 

5(13,* 

) ' 

upper : 


' , phiO 

write ( 13, * 

) ' 

lower : 


' , phim 

write ( 13, * 

) ' 

Lengths : 


1 

write 

5(13,* 

) ' 

upper : 


' , D0L 

write ( 13, * 

) ' 

lower : 


' , DmL 


end if 

write (13,*) 'Maximum nominal deviation between NURBS curve' 
write(13,*) 'and ice data: ' , epmax 

write (13,*) 'Nominal rms deviation between NURBS curve' 
write (13,*) 'and ice data: ' , rms 

close (13) 

open ( 14 , f ile= ' cukb ' / /run/ / ' . dat ' ) 

write(14,60) (CA (1, i) , CA (2, i) , i=0,mdata) 

mdata=m2datamax 

du=l . 0/dble (mdata) 

u ( 0 ) =0 . 0 

do i=l, mdata 

u ( i ) =u ( i-1 ) +du 
end do 

call FindSpanA (mdata, n, p, u ,U,spanA) 
call NCurveA( mdata, n, p, pmax, spanA, u, 

1 P, U,CA ) 

open ( 7 , f ile= ' icecv ' / /run/ / ' . dat ' ) 
do i = 0, mdata 

write (7, 10) CA (1, i) , CA (2, i) 
end do 
close ( 7 ) 

stop 

10 format (lx, 3 (e22 . 16, 2x) ) 

20 format (lx, i3, 3x, f 14 . 10) 
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30 format (lx, i3, 3x, f 14 . 10, 3x, f 14 . 10) 

40 format (lx, 10i5) 

50 f ormat ( lx, f 6 . 1 , ' and ',f6.1) 

60 format (lx, 2 (f 14 . 10, 2x) ) 
end 


integer function FindSpan (n, p, u, U) 


FindSpan calculates the span index of a curve parameter u via a 
bisection routine 


Record of Revisions: 


Date 


Programer 


Description of Change 


05/20/03 


Loren H. Dill 


Original code 


Main Variables 

n Largest index of control pnts vector 

p Degree of basis functions 

u curve parameter 

U knot vector 


Reference: Piegl & Tiller, p.68 


implicit none 

integer n, p, low, high, mid 
double precision u, U(0:n+p+l) 

if (u .eq. U(n+1)) then 
FindSpan = n 
return 
end if 

low=p; high=n+l ; mid= (low+high) /2 
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do while ( (u .lt.U(mid)) .or. (u . ge . U(mid+1))) 
if ( u .it. U (mid) ) then 
high = mid 
else 

low = mid 
endif 

mid = (low+high)/2 
end do 

FindSpan = mid 

return 

end 


^i^^i^^i^^i^i^-k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

subroutine FindSpanA (mdata, n, p, ukb, U, spanA) 

icicifificicicicici^ifififificific'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

* FindSpanA calculates the span indices corresponding to an array 

* of curve parameters ukb. Uses same algorithm as FindSpan, but 

* modified for an array of curve parameter values. 
************************************************************************ 

* Record of Revisions: 

• k 

* Date Programer Description of Change 

* 05/20/03 Loren H. Dill Original code 

■k 

************************************************************************ 


•k 
: k 
* 
* 
: k 
* 
: k 
k 
' k 
' k 
* 
* 


Main Variables 


INPUT 

mdata 

n 

P 

ukb 

U 


Largest index of data array ukb 
Largest index of control pnts vector 
Degree of basis functions 
1-D array of curve parameter values 
knot vector 


OUTPUT 

spanA 1-D array of indices for ukb data 


i k 


~k 
■k 
: k 
■k 
■k 


■k 


* Reference: Piegl & Tiller, p.68 
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implicit none 


integer i, mdata, n, p, low, high, mid 
integer spanA (0 :mdata) 

double precision ukb (0 rmdata) , U(0:n+p+l) 
do i=0, mdata 

if(ukb(i) .eq. U(n+1)) then 
spanA (i) = n 
else 

low=p ; high=n+l ; mid= ( low+high) / 2 

do while ( (ukb ( i ) .lt.U(mid)) .or. (ukb(i) . ge . U(mid+1))) 
if ( ukb(i) .It. U (mid) ) then 
high = mid 
else 

low = mid 
end if 

mid = (low+high) /2 
end do 

spanA (i) = mid 
end if 
end do 

return 

end 


k k 

subroutine NBasis ( i, mknot, p, u, U, N ) 

k k 

* NBasis calculates the p+1 non-zero basis functions within knot * 

* span index i * 

k k 

* Record of Revisions: 


Date 


05/20/03 


Programer 


Loren H. Dill 


Description of Change 


Original code 
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INPUT Variables 
i knot span index of u 

mknot largest index of U 
p degree of NURBS curve 

u value of curve parameter 

U knot vector of length m+1 

OUTPUT Variable 

N Basis function array: (N(i-p,p), N(i,p)) 


-k 


~k 

■k 


■k 
■k 
■k 
: k 
: k 


* Reference: Piegl & Tiller, p.70 

■k 


implicit none 


integer i, j , mknot, p, r 

double precision u,U(0 : mknot) ,N(0:p),left(p), right (p) , saved, temp 


N (0) =1 . 0 
do j=l,p 

left(j) = u-U(i+l-j) 
right(j) = U(i+j)-u 
saved = 0.0 
do r=0, j-1 

temp = N ( r) /( right ( r+1 ) +left ( j-r) ) 

N(r) = saved + right (r+1) * temp 

saved = left (j-r) * temp 

end do 

N ( j ) = saved 
end do 
return 
end 

■k : k 

subroutine NBasisA( mdata, mknot, p, pmax, spanA, ukb, U, NA ) 


■k'k-k'k'k'k'k-k'k'k'k'k'k'k'k'k'k'k'k'k'k'k-k'k'k'k'k'k'k'k'k'k'k'k-k'k'k'k-k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k-k'k'k'k'k'k'k'k-k'k'k'k'k'k'k'k-k'k 

* NBasisA calculates the p+1 non-zero basis functions corresponding* 

* to each element in the 1-D array of parameter values ukb * 

•jk - ■ k 




Record of Revisions: 
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Date 


Programer 


Description of Change 


k 


* 05/2003 Loren H. Dill Original code 


k 


k 

k 

k 

k 

k 

k 

k 

k 

k 


INPUT Variables 

mdata largest index of parameter array ukb 
mknot largest index of U 

p degree of NURBS curve 

pmax highest possible degree of NURBS curve 

spanA 1-D array of span indices corresponding to ukb 

ukb 1-D array of curve parameter values 

U knot vector of length mknot+1 


k 

k 

k 

k 

k 

k 

k 

k 

k 


* OUTPUT Variable * 

* NA Basis function 2-D array. Column i contains to the * 

* p+1 non-zero basis functions for given ukb(i) * 

•jk - k 




* Reference: Piegl & Tiller, p.70 

k 


implicit none 


integer i, j , mdata, mknot !, offset 
integer p,pmax,r, spanA (0 :mdata) 

double precision ukb (0 rmdata) , U (0 :mknot) , NA (0 :pmax, 0 :mdata) 
double precision left (p) , right (p) , saved, temp 


do i=0, mdata 
NA (0, i) =1 . 0 
do j=l,p 

left(j) = ukb (i) -U (spanA (i) +l-j ) 
right(j) = U (spanA (i) +j ) -ukb (i) 
saved = 0.0 
do r=0, j-1 

temp = NA ( r, i )/( right ( r+1 ) +left ( j-r) ) 
NA(r,i) = saved + right (r+1) * temp 

saved = left (j-r) * temp 

end do 

NA ( j , i ) = saved 
end do 
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end do 
return 
end 


subroutine Cparam ( mdata,Q, ukb ) 


* Computes a normalized parametric variable based upon arc length * 

* for a 2D Cartesian curve. * 

k k 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

* Record of Revisions: 


Date 


* 4/14/03 


Programer 


Loren H. Dill 


Description of Change 
Original code 


kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

* INPUT VARIABLES 

* mdata mdata + 1 = Number of points defining curve 

* Q (x, y) Cartesian coordinates for curve (2-D array) 

k 

* OUTPUT VARIABLE 

* ukb Parametric variable for curve (1-D array) 




implicit none 
integer mdata, k 

double precision ukb (0 rmdata) , Q (2, 0 rmdata) 

Compute the parametric variable ukb based on arc length and 
normalize to unity 

ukb ( 0 ) =0.0 
do k = 1, mdata 

ukb(k) = ukb(k-l) + 


1 

sqrt ( ( 

Q (1, k) 

- Q ( 1 , k-1 ) 

)*( Q(l,k) " 

Q ( 1 , k-1 

2 

+ ( Q (2, 

k) - Q 

N) 

i 

i— 1 
* 

Q (2, k) - Q (2, 

k-1) ) 


end do 
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do k = 1, mdata - 1 

ukb(k) = ukb (k) /ukb (mdata) 
end do 

ukb (mdata) = 1.0 

return 

end 


kkkkkkkkkkkkkkkkkkk-kkk-kkk-k'k-k-k-k-k-k'k-k'k-k-k-kk-k-k-kkkk-kk-kk-k-k-k-kk-k-k-k-k-kk'kkk-k'kkkk'k-kk'k 

subroutine Knotvec ( mdata, n, p, ukb, U ) 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk'k'k'k'k'k'k'k'k'k'k'k 

* Knotvec calculates a knot vector for global approximation based * 

* upon equations (9.68) and (9.69) of Piegl and Tiller (1997) . * 

* The routine is said to guarantee that every knot span contains * 

* at least one ukb point, resulting in a matrix NTNB that is 

* positive definite and well-conditioned. * 


■k 


Record of Revisions: 


Date 


05/20/03 


Programer 


Loren H. Dill 


Description of Change 
Original code 




Main variables 


INPUT * 
mdata m+1 is number of data points * 
n n+1 is number of control points in approximation * 
p degree of nurbs curve approximation * 
ukb array of parameter values corresponding to data points * 


OUTPUT 

U knot vector of length n + p + 2 




implicit none 
integer i, j , mdata, n, p 

double precision d, alpha, ukb ( 0 :mdata) , U ( 0 : n+p+1 ) 
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do i = 0,p 
U(i)= 0.0 
U (n+p+l-i ) = 1.0 
end do 


d = dble ( mdata+1 )/dble( n-p+1 ) 
do j = 1, n-p 

i = int ( j * d ) 

alpha = dble ( j ) * d -dble ( i ) 

U(p+j) = (1.0- alpha )* ukb ( i - 1 ) + 

1 alpha * ukb ( i ) 

end do 
return 
end 

subroutine NCurve ( n,p,u, P, U,C ) 

: k ■ k 

* Subroutine NCurve calculates the (x,y) point corresponding to * 

* scalar parameter u of a nonrational NURBS 2-D curve. * 

k k 

* Record of Revisions: 

k 


Date 


Programer 


* 05/20/03 Loren H. Dill 

k 


Description of Change 
Original code 




* INPUT VARIABLES 

k 


n 


n+l is number of control points 
degree of B-spline 


* u curve parameter 

* P control points, 2-D array in column format 

* U knot vector 


•jk - k 

* OUTPUT * 
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c 


calculated (x,y) point on NURBS curve, 1-D array 


~ k -k 

implicit none 

integer i,j,k, mknot, n, p 
integer FindSpan 

double precision u, P(2,0:n ), U( 0:n+p+l ) 
double precision N(0:p) 
double precision C(2) 

mknot=n+p+l 

i = FindSpan ( n,p, u, U ) 

call NBasis ( i, mknot, p, u, U, N) 

C(l) =0.0 
C (2 ) =0.0 
j = i - p 

* Sum the nonzero terms only 

do k = 0 , p 

C ( 1 ) = C ( 1 ) + N( k ) * P( 1, j + k ) 

C (2) = C (2) + N ( k ) * P ( 2, j + k) 

end do 

return 

end 

■k 

subroutine NCurveA( mdata, n, p, pmax, spanA, u, P, U,CA ) 

: k 

■k ■ k 

* Subroutine NCurveA calculates all (x,y) points corresponding to * 

* the ID array u for a nonrational NURBS 2-D curve. 

■k ■ k 

* Record of Revisions: 

: k 

* Date Programer Description of Change 
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05/20/03 


Loren H. Dill 


Original code 


' k 


INPUT VARIABLES 


mdata largest index of parameter vector u 

n largest index of control points 

p degree of NURBS curve 

pmax maximum degree of NURBS curve 

spanA 1-D array of span indices corresponding to ukb 
u curve parameter, 1-D array 

P control points, 2-D array in column format 

U knot vector 




OUTPUT 


CA 


calculated array of (x,y) points on NURBS curve, 
2-D array 




implicit none 

integer i,j,k, mdata, mknot, n, p, pmax 
integer spanA (0 :mdata) 

double precision u(0:mdata), P(2,0:n ), U( 0:n+p+l ) 
double precision NA (0 :pmax, 0 : mdata) 
double precision CA (2, 0 :mdata) 

call NBasisA( mdata, mknot, p, pmax, spanA, u, U, NA) 
do i=0, mdata 

CA ( 1 , i ) = 0.0 
CA ( 2 , i ) = 0.0 
j = spanA (i) - p 

do k = 0 , p 


CA ( 1 , i ) 

= CA i 

(l,i) 

+ NA 

(k, i ) 

* P ( 

1, 

j + k 

CA ( 2 , i ) 

= CA i 

(2, i) 

+ NA 

(k, i ) 

* P ( 

2, 

j + k 


end do 
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end do 


return 

end 


kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 
: k 

subroutine Rkarray (mdata, n, offset, p, pmax, spanA, Q, ukb, P, U, NA, Rk) 

■k 

************************************************************************ 

* k 

* Rkarray calculates a 2-D array that contains the m-1 values for * 

* the x- and y- components for eqn 9.63 in Piegl & Tiller if end * 

* derivatives are free (offset = 1) . If end derivatives are fixed * 

* (offset = 2), the algorithm is appropriately modified. * 

************************************************************************ 

* Record of Revisions: 


* Date 


Programer 
* 05/20/03 Loren H. Dill 


Description of Change 
Original code 


k k 


* INPUT VARIABLES 

* 


~ k 
~ k 


* 
* 
* 
: k 
* 
* 
: k 
* 
: k 
k 
■k 
: k 
■k 


mdata 

n 

offset 

P 

pmax 

spanA 

Q 

ukb 

P 

U 

NA 


largest index of parameter vector u 
largest index of control points 

indicates whether endpoint derivatives are free 
(1) or fixed (2) 
degree of NURBS curve 
maximum degree of NURBS curve 

1-D array of span indices corresponding to ukb 
(x,y) prescribed data, 2-D array 
curve parameter, 1-D array 

control points, 2-D array in column format 
knot vector 

array of nonzero basis functions corresponding to ukb 


■k 


■k 


■k 


■k 
■k 
: k 


~k 




k 

* OUTPUT 


: k 
' k 
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Rk 2-D array containing m-1 values. Required for calculation 

of right-hand side of matrix equation. 



implicit none 

integer k,mdata, n, offset, p, pmax, spanA (0 :mdata) 

double precision Q (2 , 0 rmdata) , ukb ( 0 :mdata) , U ( 0 : n+p+1 ) 
double precision NA (0 :pmax, 0 :mdata) , Rk (2, 0 :mdata) 
double precision P(2,0:n) 


Calculate the Rk vectors, 
do k=l, mdata-1 
Rk ( 1 , k ) =Q ( 1 , k ) 

Rk (2, k) =Q (2, k) 
if (spanA(k) .eq. p) then 

Rk ( 1 , k) = Rk ( 1 , k) - NA(0,k) * P(1,0) 

Rk (2 , k) = Rk (2 , k) - NA(0,k) * P(2,0) 

if (offset .eq. 2) then 

Rk (1, k) =Rk (1, k) -NA (1, k) *P (1, 1) 

Rk (2, k) =Rk (2, k) -NA (1, k) *P (2, 1) 
end if 
end if 

if (offset .eq. 2 .and. spanA(k) .eq. p+1) then 
Rk (1, k) =Rk (1, k) -NA (0, k) *P (1, 1) 

Rk (2, k) =Rk (2, k) -NA (0, k) *P (2, 1) 
end if 

if ( offset .eq. 2 .and. spanA(k) .eq. n-1) then 
Rk (1, k) =Rk (1, k) -NA (p, k) *P (1, n-1) 

Rk (2, k) =Rk (2, k) -NA (p, k) *P (2, n-1) 
end if 

if ( spanA(k) .eq. n ) then 
if (offset .eq. 2) then 

Rk (1, k) =Rk (1, k) -NA (p-1, k) *P (1, n-1) 

Rk (2, k) =Rk (2, k) -NA (p-1, k) *P (2, n-1) 
end if 

Rk ( 1 , k) = Rk ( 1 , k) - NA (p, k) *P (1, n) 

Rk (2 , k) = Rk (2 , k) - NA (p, k) *P (2, n) 
end if 
end do 
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return 

end 


kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

subroutine RightHandSide (mdata, n, nmax, offset, p, pmax, spanA, NA, Rk, R) 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 
k k 

* RightHandSide calculates the R array, a n-1 X 2 array, which * 

* corresponds to the right side of eqn. 9.65 of Piegl & Tiller if * 

* end derivatives are free (offset = 1) . If end derivatives are * 

* fixed (offset = 2), the algorithm is appropriately modified. * 

k k 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

* Record of Revisions: 


* Date 

* 05/20/03 


Programer 
Loren H. Dill 


Description of Change 
Original code 


k k 


INPUT VARIABLES 


mdata largest index of parameter vector u 

n largest index of control points 

nmax largest value of n permitted 

offset flag indicating whether endpoint derivatives are free 
(1) or fixed (2) 
p degree of NURBS curve 

pmax maximum degree of NURBS curve 

spanA 1-D array of span indices corresponding to ukb 
NA array of nonzero basis functions corresponding to ukb 

Rk array corresponding to eqn. 9.65 of Piegl & Tiller 


k k 

* OUTPUT * 


Right-hand side of linear matrix equation 




implicit none 
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integer j,k,mdata, n, nmax, of £ set, p, pmax, row, spanA (0 :mdata) 
double precision NA ( 0 : pmax, 0 :mdata) , Rk (2 , 0 :mdata) , R (nmax-1 , 2 ) 
do j=of f set, n-of f set 
R(j+l-offset, 1 ) =0 . 0 
R(j+l-offset,2)=0.0 
do k=l,mdata-l 



row = j- (spanA (k)- 

-P) 



if ( (row .ge. 0) 

. and . 

(row . le. p) ) then 


R ( j + l-of f set, 1 

-n 

& 

II 

+l-of f set , 1 ) + 

1 

NA (row, k) 

* Rk ( 

M 

\ — I 


R ( j + l-of f set, 2 

) = R ( 

j+l-of f set, 2 ) + 

1 

NA (row, k) 

* Rk ( 

2, k) 


end if 




end do 




end do 
return 
end 


kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

subroutine ABMatrix (mdata, n, of f set, p, pmax, spanA, NA, NTNB) 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 
: k k 

* AMatrix computes the NTN matrix of Piegl and Tiller, p. 411, and * 

* stores the matrix in upper band form suitable for Lapack * 

* subroutine dpbsv. Matrox NTN has p superdiagonals. * 

k k 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

* Record of Revisions: 

k 

* Date Programer Description of Change 

* 05/20/03 Loren H. Dill Original code 

k 




k k 

* INPUT VARIABLES * 


k 

k 

k 

k 

k 

k 


mdata 

n 

offset 

P 

pmax 


largest index of input data 
largest index of control vectors 

flag for derivatives specification: 1 if not, 2 if 

specified 

current degree of NURBS curve 
parameter, max degree of NURBS curve 
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* spanA 1-array of span indices corresponding input data * 

* NA a 2-D array containing the p+1 non-zero basis * 

* functions corresponding to each data point. * 

k k 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

* OUTPUT * 


* NTNB The NTN (Transpose (N) * N) array using band storage * 

k k 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 


implicit none 

integer i, j , k, mdata, n, of f set , p, pmax 
integer kd, row, rowi, rowj , spanA (0 :mdata) 
double precision NA (0 :pmax, 0 : mdata) 
double precision NTNB (pmax+1 , n-1 ) 

kd=p 

do i=of f set, n-of f set 
do j=i,n-offset 

if ( max (1, j+l-of f set-kd) .le. i+l-off set ) then 
row=kd+l+i- j 

NTNB (row, j+l-offset)=0.0 
do k=l,mdata-l 


rowi = i- (spanA (k) -p) 
rowj = j- (spanA (k) -p) 

if ( (rowi . ge . 0 ) .and. (rowi .le. p) 


1 

.and. (rowj . ge . 0 

) .and. 

(rowj .le. p) ) then 


NTNB (row, j+l-offset) 

= NTNB i 

(row, j+l-offset) + 

1 

NA (rowj , k) * NA 

(rowi , k ) 



end if 




end do 
end if 
end do 
end do 


return 

end 


subroutine SpanTest (ep, mdata, n, p, Q, ukb, P , U, non span, i ) 
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SpanTest returns in variable nonspan the indices of spans in U * 
that have not converged. In order for a span j to converge, the * 
distances C (ukb (k) ) -Q (k) | must be less than or equal to ep, the 
distance criterion, for all ukb(k) located in span j. Variable i * 
gives the the number of non-converging spans in knot vector U. * 
Here, Q(k)= ( x(k),y(k) ) is one of the prescribed iced-airfoil * 

data points. * 


■k 


Record of Revisions: 


Date 


05/20/03 


Programer 


Loren H. Dill 


Description of Change 
Original code 




* INPUT VARIABLES * 

k 

* ep nominal tolerance between NURBS curve and each data 

* point 

* mdata largest index of input data * 

* n largest (current) index of control vectors * 

* p current degree of NURBS curve * 

* Q (x,y) the prescribed data, 2D array 

* ukb parameter array corresponding to data 

* P control points of NURBS curve 

* U knot vector of NURBS curve 

•k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 


OUTPUT 

nonspan 

i 


1-D array containing indices of non-converging spans 
number of non-converging spans 




implicit none 

integer p, n, mdata, mknot, nonspan (n+l-p) 

integer i, j , k, spanA (0 : mdata) , of f set 

double precision P (2 , 0 : n) , U ( 0 : n+p+1 ) , ukb ( 0 :mdata) 
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double precision Q (2, 0 rmdata) 

double precision ep, C(2), epsq, distancesq, dif f (2 ) 

mknot = n+p+1 
offset = 1 

epsq=ep*ep ! Actually compare the squared distances 

call FindSpanA (mdata, n, p, ukb, U, spanA) 

write ( 8 , * ) 

k=l 

i=0 

do while (k . le.mdata-1) 

call NCurve ( n,p,ukb(k), P, U,C ) 
dif f ( 1 ) = C ( 1 ) -Q ( 1 , k ) 
diff (2) = C (2) -Q (2, k) 

distancesq=dif f (1) *diff (1) +diff (2) *diff (2) 
if (distancesq . gt . epsq) then 

write (8,*) 'non-converging span, spanA(k) = ',spanA(k) 
i = i + 1 

nonspan (i)= spanA (k) 
if(spanA(k) .eq. n) then 
k= mdata 
else 

j=k+l 

do while (spanA(k) .eq. spanA(j) ) 
j = j + 1 
end do 
k = j-1 
end if 
end if 
k=k+l 
end do 
return 
end 


*********************************************************************** 

subroutine Deviations (mdata, n, p, pmax, Q, ukb, P, U, epmax, rms, CA) 
*********************************************************************** 
i k 

* Deviations returns in variable epmax the maximum nominal distance 

* between the nurbs curve and the discrete ice data by evaluating 

* all distances | C (ukb (k) ) -Q (k) | Here, Q(k)= ( x(k),y(k) ) is one 

* of the prescribed iced-airfoil data points. The routine also 
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returns curve values CA = C(ukb) for each parameter value ukb. 
Deviations also calculates and returns the root mean squared 
distance (rms) from the data to the curve. This is again a 
nominal value in the sense that the curve passes somewhat closer 
to the each data point than is represented by C (ukb (k) ) -Q (k) |. 


* Record of Revisions: 


* Date 

* 5/20/03 

* 5/22/03 


Programer 

Loren H. Dill 
L.H. Dill 


Description of Change 

Original code 
Added rms capability 


kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

* INPUT VARIABLES * 

k 

* mdata largest index of input data * 

* n largest (current) index of control vectors * 

* p current degree of NURBS curve * 

* pmax parameter, max degree of NURBS curve * 

* Q (x,y), the prescribed data, 2D array 

* ukb parameter array corresponding to data 

* P control points of NURBS curve 

* U knot vector of NURBS curve 




* OUTPUT 

* 


~ k 
k 


k 

k 

k 

k 

k 

k 

k 


epmax 

rms 


CA 


maximum nominal deviation between NURBS curve and 
prescribed data 

root mean square deviation between the curve and 
prescribed data 

array containing all the x-y values along the NURBS 
curve corresponding to curve parameter values ukb 

k 




implicit none 

integer p, n, mdata, mknot , pmax 
integer k, spanA (0 rmdata) 

double precision P (2 , 0 : n) , U ( 0 : n+p+1 ) , ukb ( 0 :mdata) 
double precision Q (2,0 rmdata) 

double precision epmax, epmax2, rms, sum, CA (2, 0 rmdata) 
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double precision dist2, diff (2) 


mknot = n+p+1 
epmax2=0 . 0 
sum=0 . 0 

call FindSpanA (mdata, n, p, ukb, U, spanA) 

call NCurveA( mdata, n, p, pmax, spanA, ukb, P, U,CA ) 

do k=l,mdata-l 

diff (1)= CA (1, k) -Q (1, k) 

diff (2)= CA (2, k) -Q (2, k) 

dist2=diff (1) *diff (1) +diff (2) *diff (2) 

sum=sum+dist2 

epmax2=dmaxl (epmax2, dist2) 
end do 

epmax = sqrt (epmax2 ) 

rms = sqrt ( sum/dble (mdata-1 ) ) 

return 

end 


kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

subroutine NewU (mdata, n, nonspan, inonspan, p, ukb, U, Unew) 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 
• k 

* NewU adds a knot at the midpoint of every span that contains 

* data points that have not met the distance criterion. 

* 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

* Record of Revisions: 

: k 

* Date Programer Description of Change 

* 05/20/03 Loren H. Dill Original code 

i k 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 


■k 

k 

* INPUT VARIABLES 

k 

* mdata Largest index of ukb, the parameter array 

* n On input, largest index of sum in current NURBS curve 

* nonspan 1-D array containing the span indexes of current knot 

* vector that have not converged. 
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* inonspan number of elements in nonspan 

* p degree of NURBS curve 

* ukb the parameter array, of length mdata + 1 

* U Original knot vector 

* OUTPUT VARIABLE 

i k 

* n Largest index of sum in new NURBS curve 

* Unew New knot vector 

' k 

implicit none 

integer i, inonspan, j , mdata, mknot, n, nonspan ( inonspan+1 ), p 

double precision ukb (0 :mdata) , U (0 :n+p+l) 

double precision Unew (0 : n+p+l+inonspan) , unew (inonspan) 

* Set nonspan ( inonspan+1 ) to mknot for termination 

* condition 

mknot=n+p+l 

nonspan (inonspan+1) =mknot 

do i=0,p 

Unew ( i ) =0 . 0 

Unew (mknot + inonspan-i ) =1 . 0 
end do 

do i=l, inonspan 

unew(i) = 0.5*( U( nonspan(i) ) +U ( nonspan(i) + 1) ) 

end do 

j = 1 

do i = p+1 , mknot+inonspan-p-1 

if(i .le. nonspan(j)+ j -1) then 
Unew(i)= U(i -j + 1 ) 
else 

Unew ( i ) =unew ( j ) 
j = j + l 
end if 
end do 

n = n + inonspan 
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return 

end 


************************************************************************ 

subroutine Verify (mdata, n,p,ukb,U, success) 
************************************************************************ 

* 

* Verify checks knot vector U to insure at least one parameter 

* value is located within each knot span. If it does, output 

* variable success reports True. If not, success reports False 

* 

************************************************************************ 

* Record of Revisions: 

* 

* Date Programer Description of Change 

* 05/20/03 Loren H. Dill Original code 

* 

************************************************************************ 


* INPUT VARIABLES 

* mdata Largest index of ukb, the parameter array 

* n Largest index of control points 

* p degree of NURBS curve 

* ukb the parameter array, of length mdata + 1 

* U Original knot vector 

* 

************************************************************************ 

* OUTPUT VARIABLE 

* 

* success Logical variable indicating if .true, that a element of 

* ukb is located within each span of new knot vector 
************************************************************************ 


implicit none 

integer i, mdata, n, p , spanA (0 :mdata) 
double precision ukb ( 0 :mdata) , U ( 0 : n+p+1 ) 
logical success 

success = .true. 

call FindSpanA (mdata, n, p, ukb, U, spanA) 


do i=l,mdata-l 
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if ( spanA ( i+1 ) -spanA ( i ) . gt . 1 ) then 

success = .false. 

write (8,*) 'No ukb element in span spanA (i)+l 
return 
end if 
end do 

return 

end 
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